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PREFACE 


This document contains a set of summary presentations of research in 
progress at the Langley Research Center in Computational Fluid Dynamics (CFD). 
While this document does not include all or even most of the work going on in 
CFD at Langley, it does display in summary form the general scope and nature 
of research in this area. 

Computational Fluid Dynamics is, of course, an important and substantial 
research activity at the Langley Research Center. No single organization has 
the primary responsibility for the CFD effort at the Center, but rather 
research activities in CFD are carried out within many different organiza- 
tional units as an integral part of a wide variety of research programs. 

The purpose of this compendium is to help identify, through numerous 
summary examples, the scope and general nature of the CFD effort at Langley. 
This document will help inform researchers in CFD and line management at 
Langley of the overall CFD effort. In addition to the in-house efforts and 
work at ICASE and JIAFS, out-of-house CFD work supported by Langley through 
industrial contracts and university grants are included. Researchers were 
encouraged to include summaries of work in preliminary and tentative states 
of development as well as current research approaching definitive results. 


Wayne D. Erickson 

Senior Scientist, Office of the Director 


Jerry C. South, Jr. 

Head, Theoretical Aerodynamics Branch, STAD 


Jerome T. Foughner, Jr. 

Technical Assistant to the Senior Scientist 
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1. Sub/Supersonic Linear Theory Methods 


Methods which yield solutions to aerodynamic problems by superposition of 
fundamental solutions of the Prandtl -Glauert linear equation for the potential. 



DEVELOPMENT OF A FREE VORTEX SHEET THEORY FOR ESTIMATION 
OF AERODYNAMIC LOADS ON WINGS WITH SEPARATION INDUCED VORTEX FLOWS 

Janies M. Luckring* 


Separation induced vortex-flows from the leading and side edges play an 
important role in the high angle-of-attack aerodynamic characteristics of a 
wide range of modern aircraft. In the analysis and design of high-speed air- 
craft, a detailed knowledge of this type of separation is required, particu- 
larly with regard to critical wing loads, optimized vortex lift characteristics, 
and the stability and performance at various off-design conditions. The pur- 
pose of this research is to develop an analytical method based on advanced 
panel technology capable of estimating three-dimensional pressure distributions 
for wings with separation induced vortex flows. The code development is being 
performed by the Boeing Company, Seattle, Washington, under contract to NASA 
Langley. Numerical and experimental evaluation studies are being conducted 
jointly. 

Figure 1 shows a typical panel arrangement and a brief description of the 
method follows. The wing, free sheet, and near wake are paneled with biquad- 
ratic doublet singularities. Thickness effects are modeled with bilinear sur- 
face source singularities. The fed sheet is a simplified model of the vortex 
core and is merely a kinnematic extension of the free sheet. Zero mass flux 
boundary conditions are imposed on the wing and free sheet; zero pressure jump 
boundary conditions are imposed on the near wake. Because the vortex sheet 
must be locally force free, the zero pressure jump boundary condition is addi- 
tionally imposed. As a consequence, the resulting system of equation is non- 
linear (both the strength and shape of the free vortex sheet are unknown) and 
an iterative solution procedure is necessitated. 

A fairly wide range of application of this method has been conducted. In 
figure 2, comparison of the experimental surface pressure distributions with 
the free vortex sheet theory as well as attached flow theory is presented for 
a delta wing at angle of attack and sideslip. The free vortex sheet theory 
provides a reasonable estimate of the experimental loads. Also apparent in 
figure 2 is the ability of the theory to predict the appreciable nonconi cal 
nature of the flow as well as the inadequacy of attached flow methods. 

Current research efforts are focused on resolving convergence sensitivi- 
ties encountered for additional geometries of interest such as cropped wings, 
wing-body configurations, and wings with appreciable twist and camber. 

Additional efforts are directed at improving code efficiency. 


*STAD, 505-31-43, (804) 827-3711 
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Figure 2.- Spanwise pressure distributions in sideslip 





EXTENSION OF THE FREE VORTEX FILAMENT 
THEORY TO INCLUDE CORE MODELING 


James F. Campbell 


Recently, lifting-surface methods have been applied to calculate the 
aerodynamics of wings having leading-edge vortex-flows. Two popular methods 
use doublet-panels and free- vortex filaments. The panel method can predict 
accurate results, but has not been extended to treat complex wing geometries 
which can have multiple vortex systems. Although the free-vortex filament 
method is easier to extend to the more complex geometries, its predicted 
pressure distributions are not as accurate because the free vortices are more 
diffused than the concentrated vortex found experimentally. 

The current effort is being conducted by Dr. E. C. Lan of the University 
of Kansas to improve the pressure predicting capability of Mehrotra's free 
vortex filament model (ref. 1) by allowing the filaments to merge into a 
concentrated vortex core. Mehrotra's model is used to start the solution, 
then the centroid of the vortex filament system determines the vortex core. 
The boundary conditions are: (1) no flow through the wing, (2) Kutta 

conditions at leading and trailing edges, and (3) force-free conditions on 
the free vortex filaments, concentrated core, and trailing wake elements. 
Figure 1 shows a sketch of the vortex filament model and figure 2 illustrates 
the improved pressure predictions after the core is added. The core model 
also obtains the converged solution in about one half the computer time. 

The method is currently being exercised on cambered wings and being 
extended to the wing-strake arrangement having multiple vortex cores. 

Reference: 

1. Mehrotra, S. C.; and Lan, C. E.: "A Theoretical Investigation of the 

Aerodynamics of Low-AspectrRatio Wings With Partial Leading-Edge 
Separation," NASA CR-145304, January 1978. 
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Figure 1.- Sketch of Mehrotra's vortex filament model. 


Theory 

Vortex Filament Model (Mehrotra, 1978) 



Figure 2.- Effect of leading-edge vortex core j model on pressure predictions 
for a delta wing. 
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THE DEVELOPMENT OF SUPERIOR SINGULARITIES FOR MODELING POTENTIAL 
FLOW AERODYNAMIC PROBLEMS AT SUBSONIC AND SUPERSONIC SPEEDS 

F.A. Woodward^ 2 
Charles H. Fox, Jr.f 
Emma Jean Landrum"^ 

The potential flow aerodynamic programs available in the mid-1960's 
employed relatively simple singularity models. In 1969, analytical studies 
were initiated to support the development of improved singularity models for 
use in a unified subsonic-supersonic potential flow aerodynamic program. This 
computer program, USSAERO, has been used as the test bed for these analytical 
studies. The USSAERO program has also been widely distributed in the United 
States as well as in Europe and South Africa. The currently released versions 
of USSAERO utilize various combinations of constant and linearly varying 
source and vortex distributions. 

After 10 years of working with source and vortex singularities, it was 
discovered that a very special combination of them resulted in something with 
highly desirable properties. This significant discovery is considered a new 
singularity concept and is called a triplet. The triplet possesses numerous 
analytical properties which make it superior to other singularities. These 
properties are currently being explored. A computer program utilizing 
triplet singularities should have significant advantages over existing 
programs. 

References: 

1. Contract, NASl-10408 "Analytical Study of Aerodynamic Characteristics of 
Wing-Body-Tail Combination at Subsonic and Supersonic Speeds." 

2. Contract, NASl-12900 "Calculation of Aerodynamic Characteristics of Wing- 
Body-Tail Configurations in Subsonic, Transonic, and Supersonic Flow." 

3. Contract NASl-15792 "Generalized Triplet Singularity, Analytical Methods, 
Inc." 

4. Woodward, F. A.: "An Improved Method for the Aerodynamic Analysis of Wing 

Body-Tail Configurations in Subsonic and Supersonic Flow, Part II - Com- 
puter Description," NASA CR-2228, Part II, May 1973. 

5. Woodward, F, A.: "USSAERO Computer Program Development, Versions B and C, 

April 1980, NASA CR-3227. 

6. Winter, 0. A.: "The Incorporation of Plotting Capability into the Unified 

Subsonic Supersonic Aerodynamic Analysis Program, Version B," NASA 
CR-3228, 1980. 

7. Fox, C. H., JR and Breedlove, W. J., Jr.: "Application of an Improved 

Unified Subsonic-Supersonic Potential Flow Method for the Aerodynamic 

Analysis of Aircraft Configurations," AIAA Paper No. 74-186, 1974. 

8. Woodward, F. A. and Landrum, E. J.: "The Supersonic Triplet-A New Aero- 

dynamic Panel Singularity with Directional Properties," AIAA Paper No. 79- 
0273R, Vol . 18, No. 2, Page 138, February 1980. 

9. Landrum, E. J. and Miller, D. S.: "Assessment of Analytic Methods for the 

Prediction of Aerodynamic Characteristics of Arbitrary Bodies at Super- 
sonic Speeds," AIAA Paper No. 80-0071, 1980. 
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EXTENSION OF TWO-DIMENSIONAL ANALYSIS AND DESIGN METHOD 
TO THREE DIMENSIONS 


James L. Thomas* 

At subsonic speeds, the only general inviscid calculation methods 
presently available are surface panel methods which distribute source and/or 
doublet singularities on the surface geometry. An analysis and design 
formulation, which has proven successful in two dimensions, is being extended 
to allow three-dimensional design and analysis applications. This work is 
being performed under contract to NASA Langley by the McDonnell Douglas 
Corporation. 

The two-dimensional analysis method uses linear distributions of source 
and vortex singularities on planar panels with normal velocity boundary 
condition satisfaction achieved by an internal perturbation potential boundary 
condition. The method is numerically stable and the prediction accuracy for 
surface velocity calculated from singularity strength is competitive with more 
complex curved panel formulations. 

The design method is iterative and uses the derivative matrix of velocity 
with respect to geometry perturbations in the inversion calculations. Since 
such an approach accounts for the underlying physics of the flow, the design 
iteration converges very rapidly. Figure 1 illustrates the rapid convergence 
of the design cycle for a two-dimensional airfoil with prescribed pressure 
distributions available from an exact conformal method solution. 

Figure 2 shows an application of the analysis cind design method to an 
airfoil with trailing-edge separation. The inviscid model is based on 
semi -empirical rules for the behavior of pressure distributions in separated 
flow. The separation point and velocity level at a downstream point are spe- 
cified and the analysis and design program (MAAD) iterates for the wake shape. 
The simple model predicts the pressure distribution in good agreement with the 
experimental data. 

The three-dimensional extension will incorporate the two-dimensional 
design formulation into a three-dimensional method which uses a least square 
quadratic doublet formulation similar to the two-dimensional linear vortex 
formulation. 


»*STAd 7 5b5-3i-^^ 8011-827-3611 


8 





After 2 Iterations 



Figure 1 Iteration cycle for two-element airfoil 

inverse design. 



Figure 2.- Application of potential-flow combined analysis 
and design program to separated airfoil analysis. 



DEVELOPMENT OF A VECTOR-CONTINUOUS LOADING CONCEPT FOR 
AERODYNAMIC PANEL METHODS 


James L. Thomas* 

For computation of inviscid flow over aerodynamic configurations at sub- 
sonic and supersonic speeds, linearized aerodynamic influence methods using 
surface panel methods have been used for over a decade with reasonable success 
even for analyses of very complex shapes. Much effort has been expended on the 
development of such methods to minimize discretization errors and problems with 
solution stability. Such problems have rendered many panel method formulations 
unreliable unless used with considerable finesse. A new approach to the reduc- 
tion of discretization errors in such methods has been developed (ref. 1) and 
applied to two-dimensional incompressible analyses. The method is being 
contrasted with higher order curved panel formulations and also being studied 
for further development to three-dimensional applications. 

The approach is based on preventing the occurrence of induced velocity 
singularities at panel slope discontinuities by maintaining continuity of the 
velocity jump vector, defined as the vector sum of panel normal and tangential 
velocity differences arising from the source and vortex singularities, 
respectively. The sketches in figure 1 illustrate several features of the 
loading vector concept. Triangular loadings of source and vorticity for a 
single panel are shown in figure 1(a). Loading vector continuity shown in 
figure 2(a) indicates that discontinuous source and vortex strengths are 
required across nonparallel panels. The individual source and vorticity 
strengths on a panel can be described by a general loading vector at panel 
intersections shown in figure 1(c). The overlapping of adjacent panels is 
illustrated in figure 1(d) and shows the continuous variation in both direction 
and magnitude of the loading vector along a panel. 

The approach has been evaluated by application to several external and 
internal flow problems. Figure 2 illustrates the error trends for a circular 
cylinder analysis using the load vector continuity formulation and a conven- 
tional source, vortex continuity formulation. The loading vector formulation 
exhibits a second order accuracy trend even in the prediction velocity at panel 
corners which becomes singular in the conventional source, vorticity continuity 
formulation. In other applications, the method was shown to be generally well- 
behaved and essentially insensitive to irregularities in panel size 
distribution. 


References: 

1. Kemp, W. B. , Jr.: A Vector-Continuous Loading Concept for Aerodynamic 

Panel Methods. NASA TM-8010»I, May 1979. 
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Figure 1.- Properties of the vector-continuous loading 
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Figure 2.- Discretization error trends in circular 
cylinder analysis. 



ROTOR BODY INTERFERENCE 


Carl E. Freeman* 

The fuselage of a single main rotor helicopter operates in a flow field 
which is composed of its flight velocity, ambient wind conditions, induced 
flow from the engines, and induced velocities from the rotors. Likewise, the 
main rotor system is affected by the first three factors and induced veloci- 
ties created by the fuselage altering the free-stream conditions. 

The trend in modern helicopters is to decrease rotor-fuselage separation 
in order to increase transportability and decrease the overall vehicle size. 
Accompanying the decrease in rotor-fuselage separation has been a requirement 
to reduce the fuselage vibration level to provide a more comfortable cabin and 
cockpit, especially in nap-of-the-earth flight. The conflict is self evident. 
As separation has decreased, dynamic blade loads and rotor-induced pressure 
variations on the fuselage have increased. Therefore, the analysis of the 
diverse flow field at any condition from hover to high-speed flight is impor- 
tant in helping to minimize the dynamic blade and fuselage loadings. 

Even though the dyneimic nature of the rotor wake is important in 
analyzing certain problems such as vibration, there are many design details 
which can be analyzed using a time-averaged downwash field. Therefore, a 
computer code has been developed which combines a vortex tube rotor-wake 
theory with an incompressible, potential-flow panel method. This method 
calculates both on-body and off -body velocities, fuselage surface pressures, 
and total loads. 

Figure 1 presents flow-field velocity vectors as calculated by the 
combined vortex-tube/panel method for a low-speed flight condition. Figure 2 
presents experimental and calculated time-averaged surface pressures on the 
top center line of the fuselage. At this low-speed flight condition, correla- 
tion is good. 


<*STAD, 505-42-13, 80l|-827-36ll 
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APPLICATION OF SUBPANELING CONCEPTS TO THE MODELING OF 
VORTEX/SURFACE INTERACTIONS 

W. Elliott Schoonover, Jr.* 

Many leading-edge-vortex flow situations are difficult to model with 
panel methods because they require close proximity of the panels of the vortex 
sheet and those of the wing surface. Close proximity relative to panel size 
leads to numerical instabilities and divergence of iterative solutions. 
Examples of such situations include leading-edge vortices at low angles of 
attack, vortex separations from rounded surfaces and the secondary-flow re- 
gions between a leading-edge vortex and the wing leading-edge. Circumventing 
this difficulty by simply increasing panel density quickly becomes prohibi- 
tively expensive in terms of computational costs. A technique is being devel- 
oped which will allow localized increasing of panel density in those flow 
regions where increased density is required. This technique, called sub- 
paneling (ref. 1), is expected to economically solve many current problems in 
vortex flow modeling. 

An example of the use of subpaneling is shown in the attached figure. 

In the lower left is shown a Clark-Y airfoil at incidence with a vortex fixed 
above its upper surface by a sink. The resulting pressure distribution is 
shown in the upper left. The solid line is an exact solution developed by 
Rossow (ref. 2), using a transformation technique and is shown as a benchmark 
for evaluating the improvements made by subpaneling. The circular symbols in- 
dicate calculated pressures obtained by modeling the airfoil with 40 panels. 
The right half of the figure shows a detailed scan of the upper surface pres- 
sures in the vicinity of the suction peak beneath the vortex. Although the 
pressures elsewhere around the airfoil were calculated quite accurately, with- 
out subpaneling (circular symbols), the suction peak is substantially under- 
estimated. The square symbols indicate a more accurate solution obtained by 
subpaneling; specifically, each panel within four panel widths of the calcula- 
tion point is replaced by five subpanels. Each subpanel control point lies 
on an interpolation of the airfoil contour and each subpanel singularity 
strength is a biquadratic interpolation of a set of surrounding panel singu- 
larity strengths. The improvement in accuracy is dramatic. The key point, 
however, is that this accuracy improvement was obtained after the forty panel 
singularity strengths had been determined and, thus, at a cost significantly 
less than by simply modeling the entire airfoil with 200 panels. 

All effort to date has been two-dimensional; the long-range objective, 
however, is extension of subpaneling techniques to the 3-dimensional paneling 
used to model leading-edge vortex flows. This research is being conducted by 
Analytical Methods, Inc., Bellevue, Washington, under NASA contract NASl-15495. 

References : 

1. Maskew, Brian: "A Subvortex Technique for the Close-Approach to a 
Discretized Vortex Sheet." NASA TM X-62487, September 1975. 

2. Rossow, V. J.: "Lift Enhancement by an Externally Trapped Vortex, Journal 
of Aircraft . Vol . 15, No. 9, September 1978, pp. 618-625. 
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Figure 1.- Improvement of Calculated Pressure Distribution Accuracy Through Subpaneling. 



INVISCID ANALYSIS OF TIP VORTEX FORMATION 
ON OSCILLATING WINGS 
Warren H. Young, Jr.* 


The lack of analytical design tools for helicopter rotor blade tips has 
hindered the development of improved tip planforms for decreased control loads, 
reduced rotor noise, and optimized performance. This program is one of the 
required tools. The purposes of this analysis are studies of competing tip 
shapes and to provide guidance for improved tip vortex modeling in future 
rotor wake programs. 

The analysis (ref. 1) is an unsteady potential flow panel method using 
planar quadrilateral panels to represent the surface of thick wing tips 
(fig. 1). Each panel has a constant source and doublet distribution. The 
Dirichlet boundary condition on velocity potential is applied at a central 
control point. The method includes detailed paneling around the tip edge. As 
the calculation is stepped in time, a wake separation line may be specified 
on the tip. Both the usual unsteady wake shed from the trailing edge and the 
wake shed from the tip are convected downstream at the local velocity. This 
gives the time history of the rollup of the wake to form the tip vortex. 

Fourier analysis is used to find the real and imaginary coefficients (fig. 2) 
after a complete angle of attack oscillation has been calculated. 

The analysis is designed to calculate rigid body pitching motion of an 
unstalled wing in incompressible flow. The resolution of the details of the 
tip vortex requires several hundred panels and 30 to 50 time steps. Another 
limitation of the program is that the position of the tip separation line 
(fig. 1) is an input. Three measures to avoid error in placing the seperation 
line are planned. First, the calculated pressures will be correlated to 
unsteady surface pressure measurements in the DFVLR 3 x 3-meter wind tunnel 
for four tip shapes (fig. 3). Second, an integral boundary layer code will be 
added to the analysis. Third, a quasi -steady parabolized Navier-Stokes solu- 
tion will be developed under a separate program. It will calculate the 
boundary layer, free shear layer, and vorticity diffusion. 

The development of this analysis is continuing at Analytical Methods, Inc. 
under contract NASl-15472. In addition to the boundary layer, plans include 
a compressible approximation by restricting the magnitude of the angle of 
attack variation, motions other than pitching, and wind tunnel wall corrections. 

References : 

1. Maskew, B. : Influence of Rotor Blade Tip Shape on Tip Vortex Shedding-- 

An Unsteady, In viscid Analysis. AHS Paper 80-6, May 1980. 
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STEADY, OSCILLATORY, AND UNSTEADY SUBSONIC AND SUPERSONIC AERODYNAMICS 

(SOUSSA) 

E. Carson Yates, Jr.* 


Under NASA grant. Dr. Luigi Mori no and his associates at Boston 
University have developed the theory and methodology for a very general 
potential -flow analysis that can generate pressure distributions and loads 
on complete aircraft or other bodies having arbitrary shapes, motions, and 
deformations. The analysis generated an integral equation for the velocity 
potential by application of Green's theorem to the governing partial 
differential equation. This integral equation has been implemented for 
computation by surface-panel discretization and Laplace transform time 
solution. The method has been incorporated into a computer program called 
SOUSSA (Steady, Oscillatory, and Unsteady Subsonic and Supersonic Aero- 
dynamics). The program is flexible, user-oriented, and modular to facilitate 
inclusion of new capabilities as they become available (e.g., transonic flow 
which is under current development). The program should be a versatile 
and effective tool for aerodynamic, aeroelastic, and stability analyses as 
well as for structural -design optimization. 

Rockwell International recently used SOUSSA aerodynamics in a 762-panel, 
30-mode subsonic flutter analysis of the entire shuttle orbiter. Both 
symmetric and antisymmetric flutter were investigated with modes of motion 
that included rudder and eleven rotation as well as rigid-body motions and 
structural deformations. The entire external surface of the vehicle was 
paneled, as shown in the attached figure, so that finite thickness of 
lifting surfaces and aerodynamic interferences were accounted for. The 
figure alsoshows some representative results for antisymmetric flutter which 
involved primarily vertical-tail motion. Comparison with the flutter 
boundary calculated with aerodynamics for the isolated vertical tail shows 
that the aerodynamic interference of the orbiter on vertical -tail flutter 
is favorable. The orbiter is indicated to be flutter-free by a substantial 
margin. 


*Aeroelasticity Branch, 533-01-13-07, 804-827-2661 
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CALCULATION OF LIFTING PRESSURES AND 
GENERALIZED FORCES ON A THIN LIFTING SURFACE 
WITH LEADING- AND TRAI LING-EDGE FLAP-TYPE 
CONTROLS UNDERGOING HARMONIC MOTIONS 

H. J. Cunningham* 


Analysis is made, based on linear potential-flow theory, by a subsonic- 
kernal-function method with downwash collocation on the surface. Each term 
in the lifting pressure series is smooth and continuous on the surface 
except for appropriate singularities at edges and hinge lines. 

Leading and trailing edges can each be made up of from one to nine 
straight-line segments. The edges are continuous, but can have discontinuous 
slopes, at segment juncture points. Spanwise symmetry or antisymmetry of 
motion is user specified. From zero to six controls can be accommodated. 

The user supplies the data to describe the modal deformation and motion. 
Lifting pressure differences, section generalized forces, and total 
generalized forces are calculated at user-specified points and span stations. 
The calculated results are written on formatted files that can be saved and 
subsequently input to the user's post-processor program, such as, for flutter 
analysis or plotting. 

The figure show the theoretical and experimental in-phase and out-of- 
phase lifting-pressure distributions over a wing with two trailing-edge 
controls. THe flow Mach number is M = 0, and the outer control is 
oscillating with a reduced frequency k = 0 . 372 . 
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0,0^ A experiment (AVA FB-7025) 



In-Phase Pressure Distributions Resulting from Motions of Outer 
Map. A, , = 0% = 0.66% A^ = 0% k = 0.372, and M - 0 



Out-of-Phase Chordwise Pressure Distributions Resulting from Motions 
of Outer Flap. f = 0°, ^ = 0.66°, A^= 0°, k = 0.372, and MsO 
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HYBRID VORTEX METHOD FOR LIFTING SURFACES 
WITH FREE-VORTEX FLOW 

Osama A. Kandil and Li-Chuan Chu* 

E. Carson Yates, Jr.** 


A Nonlinear Hybrid Vortex method (NHV-method) has been developed for 
predicting the aerodynamic characteristics of wings exhibiting leading- and 
side-edge separations. This method alleviates the drawbacks of the Nonlinear 
Discrete Vortex method (NDV-method, also known as the multiple line vortex 
method.) The NHV-method combines continuous-vorticity and vortex-line 
representations of the wing and its separated free shear layers. Continuous 
vorticity is used in the near-field calculations, while discrete vortex-lines 
are used in the far-field calculations. The wing and its free shear layers 
are divided into quadrilateral. vortex panels having second-order vorticity 
distributions. The aerodynamic boundary conditions and continuity of the 
vorticity distributions are satisfied at certain nodal points on the vortex 
panels. An iterative technique is used to satisfy these conditions in 
order to obtain the vorticity distribution and the wake shape. Distributed 
and total aerodynamic loads are then calculated. The method will be 
extended to unsteady flow for use in flutter analyses of wings and aircraft 
at high angle of attack. 

The attached figure illustrate the kind of vortex system that is 
modeled by the hybrid method for a simple delta wing. The lower left shows 
the shape of the right-hand side of the leading-edge vortex system at the 
longitudinal position of the trailing edge. The lower right shows flow at 
a cross section 8 percent of root chord downstream from the trailing edge. 

At this location the portion of the vortex sheet shed from the trailing 
edge has started to deform (roll up) in a clockwise sense which is opposite 
to that for a trai ling-edge vortex sheet generated in the absence of leading- 
edge flow separation. The arrows represent velocities measured by Hummel, 
and the dash-dot line indicates his estimate of the vortex-sheet location. 
Agreement with the results calculated for the low-order (coarse) paneling 
arrangement used for this illustration is good. 


*01 d Dominion University 

**Aeroelasticity Branch, 505-33-53, 804-827-2661 
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A NONLINEAR CONTINUOUS-VORTICITY MODEL 
Ali H. Nayfeh and Dean T. Mook 

Because the lifting surfaces of many modern vehicles are designed for 
supersonic flight, they generally have small aspect ratios, are highly swept, 
and have sharp leading edges and/or tips. Consequently, one must expect 
aerodynamically influential vortex systems to be generated by leading-edge 
and tip separation (these are called the leading-edge and tip wakes) even at 
moderate (4°) angles of attack or deflection. 

For a single wing without interactions, the need to model the leading- 
edge and tip wakes accurately is clearly illustrated by comparisons of the 
aerodynamic loads predicted by various lifting-surface methods with experi- 
mental data. For multi -component interactions, the need to model all the 
wakes accurately is compounded by the possibility of one component being in 
the wake of another. It is not reasonable to suppose that the strength, 
direction and position of the wake can be assigned accurately by employing 
some general assumptions; rather, one should expect the solution to provide 
these. 

In earlier papers, we developed an accurate, discrete- vortex model of 
the wakes and coupled this with a vortex-lattice model of the lifting surface. 
This compound model, which provides the strength, direction and position of 
the vorticity in the wake, can be used to predict accurately the aerodynamic 
loads on far-coupled and close-coupled wings, wing-body combinations and 
wings in unsteady flows. 

Though this discrete-vortex technique yields the total loads accurately, 
it does not always lead to accurate predictions of the pressure. Thus, a 
continuous-vortici ty technique is being developed to treat lifting and non- 
lifting subsonic flows. Such an approach will result in significant 
improvements in the accuracy of the predicted pressure. Then this tech- 
nique will be extended to treat general unsteady lifting flows and finally 
to modify the general unsteady technique to treat small, harmonic oscilla- 
tions efficiently. 


*Engineering Science and Mechanics Department 
Virginia Polytechnic Institute and State University 
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2. Boundary Layers 


Methods for solving the 2-D or 3-D parabolic-type boundary-layer equations. 
(Application of established boundary-layer codes to a problem would usually 
be omitted from this compendium.) 
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ADVANCED 2-D TURBULENT BOUNDARY-LAYER 


SEPARATION PREDICTION METHOD 
Suresh H. Goradia* 

Harry L. Morgan Jr.** 


For the past several years, the Lockheed-Georgia Company has been develop- 
ing methods for the prediction of separation and stall on single-component air- 
foils. For leading-edge stall, Goradia's weighted multi parameter method for 
the prediction of laminar separation and reattachment has been shown to give 
accurate stall angle-of-attack predictions for a wide variety of airfoils. 
Application of these leading-edge stall airfoils to practical aircraft wings 
have generally produced undesirable abrupt stall characteristics. The trailing- 
edge stall airfoils which are characterized by large regions of separated flow 
near the trail ing-edge have very nonlinear lift versus angle-of-attack curves 
and produce more desirable gradual wing stall characteristics. Several theo- 
retical methods are currently being developed to model this trailing-edge 
separation region by determining the effective displacement of the fluid 
streamlines in inviscid flow. These methods require an accurate prediction of 
point of turbulent boundary layer separation. 

A recent contractual effort was initiated with Lockheed-Georgia to develop 
a new method based on energy dissipation principals and the use of inviscid 
pressure distributions to predict turbulent boundary layer growth, the loca- 
tion of separation, and the viscous pressure at the point of separation. 

Using moment of momentum and energy dissipation integral equations, trans- 
formations were performed so that terms containing gradients of pressure at 
the outer edge of the boundary layer did not appear in the governing equations. 
This resulted in a set of integral equations whose solution is not jeopardized 
by large pressure gradient errors that are induced by relatively small inaccu- 
racies in the computed pressure distribution. These integral equations also 
include generalized terms which account for the effects of large velocity 
fluctuations near the point of separation and were formulated based on a set 
of laser velocimeter data obtained by Lockheed on a Wortmann-type airfoil. 
Parametric criteria for separation were then derived from dimensional analysis 
of the resulting integral equations. Several parameters were found effective 
in predicting the separation point and a large number of correlations were 
made to determine the most effective parameter. Additional parametric criteria 
were also found to accurately predict the viscous pressure at separation. The 
results of this method applied to a recently developed 17-percent thick, 
medium-speed general aviation airfoil are presented in figure 1 and show excel- 
lent agreement for both the location of separation and the viscous pressure at 
separation. 


*Lockheed-Georgia employee 

**NASA Langley employee, STAD, 505-11-13, 804-827-4514 
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SEPARATION PREDICTION FOR M$(l)-0317 AIRFOIL USING TURBULEMT SEPARATION METHOD 



CALCULATION OF HEATING RATES ON 


THREE-DIflENSIONAL CONFIGURATIONS 
H. Harris Hamilton II* 


A method has been developed for calculating laminar heating rates on 
three-dimensional configurations using the "axisymmetric analogue" (ref. 1). 

In this approach, the three-dimensional boundary- layer equations are written 
in inviscid surface-streamline coordinates and the crossflow in the boundary 
layer is assumed small, which is the case when the streamline curvature is 
small or when the wall is highly cooled. This reduces the three-dimensional 
boundary- layer equations to a form equivalent to the usual axisymmetric 
equations. Heating rates on a configuration are calculated along streamlines 
using any existing method of solving the "equivalent" axisymmetric boundary 
layer. When several streamlines are considered, the heating over the entire 
configuration can be obtained. 

The inviscid surface streamlines, metric coefficients, and boundary- layer 
edge conditions are calculated from a three-dimensional inviscid flow-field 
solution obtained from a computer code (STEIN) developed by the Giumman 
Aerospace Corporation. Coupling the present boundary- layer method with the 
three-dimensional inviscid code means that heating rates can be calculated 
on any configuration where the flow is attached and an inviscid solution 
can be obtained. 

The heating rates along the windv^ard symmetry plane of an 80" sweep 
slab delta wing at an angle of attack of 20" are presented in figure 1. 

There is good agreement between the present theory and the experimental data. 
For the same configuration, the circumferential distribution of heating rates 
at three axial stations are presented in figure 2. Again, the agreement 
between the present theory and the experimental data is very good. Similar 
results have been obtained for other configurations such as blunt cones and 
the fuselage of the Space Shuttle Orbiter (ref. 2). 


References : 
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— Present Theory 
O Experiment 


Fig, 1 Distribution of heating in windward symmetry plane 
of 80“ sweep slab delta wing; a = 20“ ; = 9.6 

Re = 3,94 X io6/m. 
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Fig, 2 Circumferential distribution of heating on 80“ sweep 

slab delta wing; a = 20“ ; = 9.6; Re^ = 3,94 x lOvni. 
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3. Viscous/Inviscid Interaction 


Methods for interacting inviscid codes with viscous (usually boundary-layer) 
codes. 
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CALCULATION OF AEROELASTIC LOADS ON WINGS 


IN TRANSONIC FLOW 

Woodrow Whitlow, Jr. 

* 

Robert M. Bennett 


The aerodynamic loads acting on flexible wings depend on the dov/nwash 
which, in turn, depends upon the elastic deformation of the wing. Thus, the 
elastic deformation, which itself is determined by the aerodynamic loads, can 
have a significant effect on wing loading. Specifically, transonic prediction 
methods that assume rigid wings--which almost all such methods do--may predict 
loads on flexible wings that differ drastically from the true values. 

Consequently, this effort, which is to be conducted at Langley Research 
Center, is aimed at coupling an aeroelastic model with a steady transonic 
prediction method to yield the flexible wing loading. The wing which is 
currently being examined is a transport wing with a supercritical section. 

That wing is shown in Figure 1, which is taken from Reference 1. The transonic 
flow field for a rigid wing is computed using the FL022 computer program^, 
which is a finite difference solution of the full potential equation. The 
wing loading obtained from that program may then be used to determine the 
aeroelastic correction to the current wing geometry. That geometry is then 
used as the new rigid wing and the process continued until a converged 
geometry and flexible wing loading is obtained. The effects of viscosity 
may also be included by incorporating the solution of the boundary layer 
equations into the iteration process. 

The type of effects that are expected are shown in Figures 2 and 3 
which are taken from Reference 1. Figure 2 shows the spanwise lift and 
pitching moment distributions on a wing operating at the cruise point, and 
Figure 3 shows those distributions for the same wing operating at a much 
higher dynamic pressure. These figures demonstrate the need to incorporate 
aeroelastic effects into the aerodynamic predictions for flexible wings. 

References 

1. Chipman, R. ; Waters, C. ; Mackenzie, D.: Numerical Computation of 

Aeroelastical ly Corrected Transonic Loads. AIAA Paper No. 79-0766. 
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Figure 1. Planform of a transport wing with supercritical sections. 



Figure 2. Predicted spanwise distri- Figure 3. Predicted spanwise distri- 
bution of lift and pitching moments at bution of lift and pitching moments 
the transonic cruise point. at a transonic low altitude point. 
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INTERACTION OF 3-D WING BOUNDARY-LAYER CODE WITH 

3-D TRANSONIC INVISCID WING CODES 

★ 

Craig L. Streett 


Transonic flow over three-dimensional wings, especially over those con- 
figurations using supercritical airfoil sections, is known to be strongly 
affected by viscous effects even when separated regions are very small. 

Both level and spanwise distribution of wing lift is changed by the decam- 
bering effect of the boundary layer. Although viscous correction using two- 
dimensional strip boundary- layer calculations shows fair agreement with 
experiment in some cases, the strong spanwise crossflow in the cove region 
of a supercritical wing dictates that any viscous correction used must be 
three-dimensional to have hope of correcting properly for these viscous 
effects. 

Iteratively interacting inviscid and boundary- layer computations to 
produce an effective body-displacement surface and its corresponding pressure 
distribution is a technique which in principle should produce results of 
adequate accuracy. However, current three-dimensional finite-difference 
boundary- 1 ayer techniques require computation times and storage of at least 
the same order as presently used full-potential inviscid calculation tech- 
niques. Thus, the cost of the many complete boundary-layer calculations 
required for interaction is prohibitive. Present three-dimensional integral 
boundary-layer techniques, on the other hand, can reduce this cost by a factor 
of 20 or more and give results in good agreement with finite-difference codes 
even in difficult flow situations. In addition, integral methods are in 
general more robust than finite-difference methods, an important consideration 
during early iterations of viscous/inviscid interaction. 

At present, a modern three-dimensional integral boundary- layer code has 
been successfully interacted with the transonic ful 1 -potential wing code 
FLO-27. Converged interacted results for three wing stations on the super- 
critical SCW-3 wing are compared with experiment in figure 1. Agreement is 
good, with much of the decambering effect of the boundary layer in the cove 
region predicted. Body effects were not accounted for in this calculation; 
work to interact the wing-body code FLO-30 with the present boundary- layer 
method is underway, and agreement is expected to improve. 


*STAD, 505-31-13, 804-827-2627 
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CALCULATION OF TRANSONIC FLOW ABOUT 
AIRFOILS WITH TRAILING-EDGE SEPARATION 
Richard W. Barnwell* 


During the past decade, computational fluid dynamics methods have been 
developed to the point that reliable predictions can be made for attached 
viscous flow past airfoils traveling at subsonic and transonic speeds. 

However, despite some pioneering efforts in the 1960's, the development of 
methods for calculating separated flow has lagged considerably. As a result, 
extensive wind-tunnel experimentation is still required to determine airfoil 
stall characteristics and control effects although design-point characteristics 
can be predicted analytically. 

Viscous airfoil flowfields. at typical flight Reynolds numbers can 
generally be approximated as inviscid outer flowfields and viscous boundary 
layers. Therefore, the development of reliable separated-flow prediction 
methods depends, in principal, on the successful coupling of existing inviscid 
and boundary-layer techniques, which have reached a high degree of development 
during the past decade. At present, several grants and contracts are being 
supported by the Airfoil Aerodynamics RTOP to develop hybrid methods of this 
type to account for trai ling-edge separation. F. A. Dvorak of AMI is incorpor- 
ating a vortex lattice representation of the separation streamlines into the 
NASA-Langley /Lockheed multi-element vortex-lattice airfoil program (MCARF) for 
subsonic flow. This technique, which is similar to the single-element method 
described in reference 1, uses one of several integral boundary- layer methods. 
L. A. Carlson of Texas A & M is developing a finite-difference/integral 
boundary-1 ayer method for both subsonic and transonic applications (ref. 2). 

0. R. Burggraf of the Ohio State University is developing a method in which 
fi nite-diference techniques are used both in the inviscid and boundary- layer 
regions. Barnwell has recently succeeded in calculating separated-trai ling- 
edge flow about an airfoil traveling at transonic speeds. The methods used 
employed a simple finite-difference potential flow method and an integral 
boundary-layer method. Results are shown in Figure 1. 

References: 
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PREDICTION OF COMPRESSIBLE TRAILING-EDGE SEPARATION EFFECTS 

NACA 0012 AIRFOIL 


theory experiment 
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COMPUTATION OF LEADING-EDGE SEPARATION 


BUBBLES ON AIRFOILS 
Richard W. Barnwell* 


Laminar boundary-layer separation can occur in the leading edge region of 
a blunt airfoil if the Reynolds number is sufficiently low that the flow 
remains laminar front the stagnation point to the minimum pressure point. 
Downstream of this peak suction point, the laminar boundary- layer separate 
almost immediately resulting in a recirculating flow region bounded by a thin 
shear layer. This separated shear layer, which is very unstable, generally 
transitions to turbulent flow so that the lower portion of the shear layer is 
energized by turbulent mixing (fig. 1). If this mixing is sufficiently 
intense, the flow will reattach quickly and a "short" bubble results. If quick 
reattachment does not occur, then a "long" bubble results. The "short" bubble 
has a relatively small effect on airfoil performance, but the "long" bubble has 
a profound effect and generally indicates the onset of leading-edge stall. The 
change from "short" to "long" bubble (called bubble bursting) may occur 
gradually or abruptly depending on the airfoil shape. 

J. E. Carter of UTSI will develop a prediction method to analyze leading- 
edge separation bubbles. In the first phase, a recently-developed subsonic 
viscous inviscid interaction procedure presented in reference 1 will be 
combined with a current airfoil analysis code and used to study the "short" 
bubble problem. In the second phase, the main thrust will be to implement an 
improved transition model and further develop this "short" bubble analysis to 
determine its potential for predicting the occurrence of bubble bursting. 
Parametric studies will be performed to determine the sensitivity of the 
leading-edge bubble to airfoil incidence shape, Reynolds number and free-stream 
turbulence level. 

The study of 0. R. Burggraf of the Ohio State University will also 
contribute to the understanding of leading-edge separation. Part of this study 
involves an analytical study of the region around a laminar separation point. 
From this analysis, insight can be gained which should lead to improved 
numerical techniques and procedures. 
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AIRFOIL LEADING EDGE SEPARATION BUBBLE 
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PATCHED VISCOUS/INVISCID INTERACTION TECHNIQUES 
FOR PREDICTING TRANSONIC FLOW OVER 
BOATTAIL NOZZLES 

Richard 6. Wilmoth* and Lawrence E. Putnam* 


A patched viscous/inviscid interaction technique is used to solve the 
subsonic/transonic flow over boattail nozzles including the effects of bound 
ary-layer displacement, boundary-layer separation, jet plume blockage and 
jet entrainment {references 1 and 2). The technique combines (see figure 1) 
the South- Jameson full potential transonic flow solution for the inviscid 
external flow with an integral boundary solution, Presz's separation model 
and discriminating streamline analysis, a shock-capturing inviscid jet 
exhaust solution and an overlaid, parabolic-marching mixing layer solution. 
Boundary-layer effects are accounted for in the conventional manner by 
adding the displacement thickness to the body geometry. For separated flows 
the discriminating streamline shape is also included to account for the 
streamline displacement caused by the separation bubble. Jet effects are 
accounted for by combining the inviscid jet and mixing layer solutions to 
define an effective viscous plume boundary which includes both the inviscid 
plume blockage and jet entrainment contributions. The entire flow field 
is solved iteratively using a viscous-inviscid underrelaxation procedure. 

A comparison with experiments of boattail pressures predicted with this 
technique are shown in figure 2. 

References : 

1. Putnam, Lawrence E.: DONBOL: A Computer Program for Predicting 

Axisymmetric Nozzle Afterbody Pressure Distributions and Drag at 
Subsonic Speeds. NASA TM-78779, 1979. 
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M < I 


INVISCIO EXTERNAL FLOW 



Figure 1. Schematic of nozzle' afterbody 
flowfield. 




number at (b) Effect of NPR at = 0.8, 

Figure 2. Comparisons between predicted and measured boattail pressure 
distributions for cold-air jet exhaust. 
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APPROXIMATE CONVECTIVE HEATING EQUATIONS 
FOR HYPERSONIC FLOW 

E. V. Zoby, J, N, Moss, and K, Sutton* 


A design problem for most reentry spacecraft is the prediction of 
convective heating rates. As a consequence of current interest in outer- 
planet entry and advanced transportation systems for Earth entry, the space- 
craft designer now needs rapid, but reliable, heating methods which must 
be capable of also assessing the effects of such problems as arbitrary 
reactive gas compositions, complex three-dimensional and/or variable-entropy 
flows, as well as the possibility of turbulent flows over large surface areas 
on the spacecraft design. 

Based on a recent investigation,^ an approximate convective heating 
method has been developed for engineering calculations of laminar and 
turbulent heating rates. The equations are applicable to nonreacting or 
reacting gas flows for either constant or variable-entropy edge conditions. 
The variable- entropy effect on the heat transfer is approximated by defining 
the boundary- layer edge properties as the inviscid values located a distance 
from the surface equal to the boundary- layer thickness. 

The results of the engineering method are in good agreement with 
available heat-transfer data, as well as boundary- layer and viscous-shock- 
layer^ (VSL) solutions. Note the comparisons of the calculated results are 
in good agreement for Earth, Venusian, and Jovian entry conditions. 

Typical comparisons of the results of the present engineering method and 
the VSL solutions are presented in Figures 1 and 2. The convective heating 
rate results are presented about a blunt 40° cone in air at = 10 in 
Figure 1 and about a blunt 45° cone at a typical Jovian entry condition in 
Figure 2. 


References: 
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Equations for Hypersonic Flows," AIAA Paper 79-1078, June 1979. 
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Figure 2.- Comparison of convective heating rate distributions for Jovian 
atmosphere. 
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4 . 


Transonic Potential Flow 


Methods for solving the irrotational approximation to transonic flow, usually 
involving the transonic velocity potential or the stream function. This 
category may often include methods, or validation of methods, which involve 
viscous/inviscid interaction as a secondary focus. 




steady Flows 


Methods designed primarily for solution of the steady problem. 
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ANALYSIS OF TRANSONIC FLOW PROBLEMS VIA 
PARAMETRIC DIFFERENTIATION 

-k 

Woodrow Whitlow, Jr. 


One of the primary difficulties associated with predicting transonic flow 
fields is that the governing eguations are generally nonlinear. Thus, 
transonic flow equations don't readily lend themselves to formal mathematical 
analysis. To circumvent the problems caused by the nonlinearity, Rubbert and 
Landahl' used the method of parametric differentiation (MPD) to linearize the 
small disturbance transonic equation. That method yielded sets of solutions 
for flow past nonlifting airfoils. 

In this study, MPD has been applied in a manner similar to that of 
Reference 1 but extended to include lifting as well as nonlifting airfoils. 
Relaxation solutions of the linearized problem were computed on a rectangular 
grid, and predictor-corrector methods were used to obtain solutions for a 
range of airfoil thicknesses and angles of attack. A typical set of nonlifting 
pressure distributions obtained using MPD is shown in Figure 1. and a compar- 
ison of a lifting pressure distribution with experimental data^ is shown 
in Figure 2. Figure 1 summarizes a set of 21 nonlifting solutions, and 
Figure 2 'summarizes seven nonlifting and two lifting solutions. Each set of 
solutions was obtained in one application MPD. 

Current efforts, being conducted at Langley Research Center, are aimed at 
utilizing MPD to examine three-dimensional steady and two-dimensional unsteady 
transonic problems. The extension of MPD to those problems is conceptually 
straightforward. In addition to thickness and angles of attack, three- 
dimensional problems offer additional parameters such as wing sweep and twist, 
and reduced frequency is a natural parameter to use for unsteady problems. 
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Figure 1 (Above). Nonlifting pressure 
distributions on biconvex airfoils. 


Figure 2 (Right). Lifting pressure 
distributions on a six percent thick 
biconvex airfoil at one degree angle 
of attack. 
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ANALYSIS OF ITERATIVE METHODS 


FOR THE TRANSONIC POTENTIAL EQUATION 
John C. Strikwerda* 


• Several mathematical aspects of iterative methods for solving the tran- 
sonic potential equation are being examined. At present, the study is con- 
centrating on a one-dimensional model problem with a simple iteration method 
so that the effect of various iteration parameters on the behavior of the con- 
vergence process can be identified. The research will be expanded in the near 
future to include two-dimensional problems and commonly used iterative methods. 
This study should provide a better understanding of how to choose appropriate 
values of the iteration parameters for several methods currently in use. 
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DESIGN OF SHOCKLESS TRANSONIC AIRFOILS 


BY THE METHOD OF COMPLEX CHARACTERISTICS 
Josd Sanz* 


The method of complex characteristics developed by Bauer, Garabedian and 
Korn has proved to be an excellent tool for the design of shockless airfoils. 
The last version of their code (Ref. 1) includes a fundamental improvement in 
the design method which consists in mapping the hodograph domain onto a circle 
thus, allowing one to design airfoils with a given input pressure distribution. 
In collaboration with NASA-Lewis and using a simple modification of the above 
code, compressor blades were able to be designed with a gap to chord ratio of 
.89 and a turning angle of 41 degrees, as shown in Figure 1. 

However, an important problem in turbomachinery technology is the ability 
to design shockless cascades with gap to chord ratios of order one half, i.e., 
solidity of up to 2. To get this capability a new code is in an advanced state 
of development. The main new concept included is a conformal mapping of the 
hodograph domain onto an ellipse rather than onto a circle. This implies a new 
method of laying paths of integration to solve the finite difference scheme 
which approximates the differential equations modeling the flow. These new 
paths of integration are more in accord with the topology that the new mapping 
introduces and are producing a substantial reduction in computing time. 

The new code is being developed with support by the NASA-Lewis Research 
Center. We expect that the new code will also be useful in the design of thin 
isolated airfoils. 
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APPLICATION AND EVALUATION OF TRANSONIC THEORIES 
FOR DESIGNING SUPERCRITICAL FIGHTER WINGS FOR TRANSONIC MANEUVER 


Michael J. Mann* 

The recent development of three-dimensional transonic computational 
methods has opened the possibility of new design procedures for aircraft 
operating at transonic conditions. The purpose of the present study has been 
to develop a theoretical procedure for the design of transonic wings. This 
procedure (ref. 1) uses the three-dimensional transonic analysis codes in an 
iterative fashion to achieve a desired pressure distribution. The method 
applies to conditions involving a large region of supersonic flow on the wing. 
Basic principles of supersonic flow are utilized to determine the necessary 
change in airfoil geometry to produce a desired change in pressure. 

The design procedure has been used to design both a swept-back wing and 
a comparable swept-forward wing for a fighter operating at transonic maneuver- 
ing conditions. The airfoils were designed to have a large region of super- 
sonic flow which terminates with a weak shock wave. The fuselage effects were 
accounted for in both designs. The swept-forward wing configuration includes 
a canard and was designed with a transonic wing-body-canard code. The 
development of the design procedure, its application, and the associated 
experimental testing have all been accomplished at the NASA Langley Research 
Center. 

The application of the design procedure to the swept-back wing, SMF-2, 
is shown in figures 1 and 2. The Mach number is .90 and the lift coefficient 
is 0.85. In figure 1, the Jameson-Caughey FLO 27 (ref. 2) three-dimensional 
transonic computer code is correlated with experimental data on an early 
version of the SMF-2 configuration. The fuselage is represented by an 
infinite cylinder and a strip boundary layer theory was used to estimate 
viscous effects. The design procedure was then used to modify the wing shape 
in order to remove the forward part of the lambda shock and reduce the 
strength of the trail ing-edge shock. The improved pressure distribution is 
compared with the original one in figure 2. These results were computed 
with the FLO 27 code. 

References: 
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Figure K- Early Configuration of SMF-2 Wing. 
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Figure 2.- Theoretical Modification to SMF-2 Wing 



IMPROVEMENTS IN THE ARTIFICIAL DENSITY METHOD FOR 


COMPUTATIONAL TRANSONICS 

* 

Jerry C. South, Jr. 


The artificial-density method (ADM) has been shown as an efficient method 
for including necessary dissipation terms in finite-difference and finite- 
element transonic-flow calculations (ref. 1). However, the method usually 
produced overshoots at mesh points just upstream of shocks; the overshoots 
were eliminated by increasing the numerical dissipation coefficient at the 
expense of more "smearing" of the shock jump. An improved procedure has 
been developed which greatly diminishes the shock overshoots without in- 
creasing the dissipation. It was found that the difficulty stemmed from the 
density calculation. 

Space does not permit presentation of the details; however, the improved 
method uses a more "compact" difference operator. In the original method, 
even at subsonic flow points, the difference operator spread over five points 
in each direction; the latest method only requires three. Some comparisons 
of the two methods are shown in figures 1 and 2. The original method is 
identified as "node-point density" and the new method as "midsegment" or 
"midcell density." Figure 1 shows results for 1-D flow in a constant-area 
duct with a normal shock. The original method gives an overshoot upstream 
from the shock while the new "midsegment density" method gives practically 
no discernible overshoot. The method has been tested extensively in 2-D 
transonic flows, an example of which is shown in figure 2 for a nonlifting, 
NACA 0012 airfoil at M = 0.85. Here again we see that the original method 
gives an overshoot at tf?e shock, while none is evident in the "midcell" den- 
sity method. Versions of the new method are now being tested in 3-D flows 
using vectorized algorithms for the CYBER 203 computer. 

References: 
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DISTANCE. X 


Figure 1.- Resolution of normal shock wave using 
two versions of artificial density 
method. 



Figure 2.- Pressure distribution for NACA 0012 airfoil. 
(M^ = .85, a = 0, 121 x 25 mesh). 
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DEVELOPMENT OF A FULLY CONSERVATIVE, AXISYMMETRIC 
TRANSONIC FLOW CODE 

•k 

Lawrence L. Green 


Technology has been developed for analysis of transonic, ax i symmetric 
bodies by South and Jameson in reference 1. A nonconservative version of the 
computer code was developed and distributed by Keller and South in reference 2. 
The present project is to apply the artificial compressibility method presented 
by Hafez, Murman, and South, reference 3, to the RAXBOD code and develop a 
fully conservative transonic flow code for axi symmetric bodies. 

Primary modifications have been successfully completed and preliminary 
results for a transonic sphere are shown in the accompanying figure. For 
comparison, the results for the same case are shown from the nonconservative 
RAXBOD code. The conservative results show increased shock strength and dis- 
placement of the shock somewhat downstream of the nonconservative version, 
as expected. Comparisons will be made in the future with the work of Salas, 
involving solutions to the Euler equations for similar problems, described 
elsewhere in this compendium. 
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DEVELOPMENT OF VECTOR ALGORITHMS FOR COMPUTATIONAL TRANSONICS 


James D. Keller* 


The state of the art in transonic flow computations has advanced to the 
point where two-dimensional flows, including boundary-layer interaction, can be 
computed efficiently on current computers. Three-dimensional calculations, 
however, are very time consuming and expensive on conventional computers such 
as the Control Data CYBER-175. If we can make proper use of the new "super" 
computers such as the CYBER-203, we should be able to reduce the cost of three- 
dimensional transonic calculations. Unfortuntely, the algorithm which has 
been used extensively in the past (successive line over-relaxation) is not 
very well -suited to the new vector processor type of computers. This was 
particularly true of the STAR-100, which had a rather slow speed on scalar 
operati ons. 

Reference 1 gives the details of a number of different algorithms which 
are not only vector! zable, but which also have a good rate of convergence. The 
one which appears most attractive so far is called ZEBRA II. It is relatively 
easy to vectorize and has a convergence rate which is the same as SLOR. The 
ZEBRA II scheme involves very few scalar operations, the shortest vector 
operations have a vector length equal to half the number of points in a cross 
plane, and the scheme has no storage problems on a virtual memory computer. 

The ZEBRA II scheme has been implemented in a small pilot code for flow over a 
swept wing using the artificial density method (ref. 2). Computer times for 
this pilot code on various computers are given in Table 1. This pilot code can 
handle only the simplest swept-wing geometry and uses Cartesian coordinates and 
planar boundary conditions. Efforts are now underway to implement the ZEBRA II 
algorithm into production-type three-dimensional codes. An in-house effort is 
underway to put the ZEBRA II algorithm into the FLO-22 code which was developed 
by Antony Jameson at New York University. This code solves the full-potential 
equation for a swept wing using a parabolic coordinate system. Computer 
Sciences Corporation has been given a task to implement the ZEBRA II algorithm 
in the WIBCO code developed by Charles Boppe at Grumman. This code solves the 
transonic small -disturbance potential equation for a swept wing using a grid 
embedding technique. 
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Table 1.- Timing comparison for ZEBRA II algorithm on 
a 3-D transonic flow code. Approximately 
30,000 grid points, 400 iteration cycles, 
density updated every 4th cycle. 


MACHINE 

CPU TIME (SEC.) 

CYBER-175 

456 

CDC 7600 

282 

CRAY-1 without vectorization 

144 

CRAY-1 with partial vectorization 

69 

STAR-100 vectorized code 

58 

CYBER- 203 vectorized code 
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BASIC RESEARCH IN CFD AT NYU 


Jerry C. South, Jr.* 


Under NASA Grant NSG 1579, the NYU Research Team has been developing 
numerical methods for transonic flows. The team presently consists of Prof. 
Paul Garabedian (director of the Fluid Dynamics Division, Courant Institute 
of Mathematical Sciences), Prof. Antony Jameson and Drs. Frances Bauer, Octavio 
Betancourt, and Geoffrey McFadden. This effort has been a continuing one since 
1970, and the past accomplishments are well known in the field of computational 
transonics: The complex-characteristics hodograph design method for shockless 

supercritical airfoils; the inviscid, circle-plane, mixed finite-difference 
airfoil analysis code; theory for iterative schemes for mixed flows; the 3-D 
yawed and swept wing codes used throughout industry and Government; the viscous 

2- D airfoil "BGK" code; the conservative full -potential equation 2-D code; and 
the 3-D finite-volume fully conservative wing-body code. In 1976, the team was 
given the NASA Group Public Service Award. 

Most recently the NYU team has focused on two problems: The transonic wing 

design problem and the multigrid relaxation method for transonic flows. The 

3- D wing design method is based partly on the 2-D design method developed by 
Dr. McFadden in his dissertation (ref. 1). For a given input pressure or speed 
distribution, the ordinates of the wing sections are adjusted iteratively in 
such a way as to move the pressures calculated from the NYU FLO-22 swept-wing 
analysis code (ref. 2) toward the desired input values of pressure. Different 
formulas are used at subsonic and supersonic points to adjust the surface 
ordinates. Preliminary results have been obtained on coarse meshes. 

The multi grid research has led to the most efficient 2-D circle-plane 
transonic airfoil code in existence. Jameson's "MAD" scheme (ref. 3) (for 
^ultigrid with ^lternating-£i recti on relaxation) achieves convergence for 
transonic airfoils in 30 to 80 cycles, with a spectral radius of about 0.6 to 

0. 8. on fine meshes. This research is being pursued further for 3-D wing 
applications. 
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COMPUTATION OF ROTATIONAL TRANSONIC FLOW USING THE STREAM FUNCTION 


Jerry C. South, Jr. 


* 


Dr. Mohamed Hafez has recently joined the GWU faculty at LaRC and will 
specialize in computational fluid dynamics. Recently, he developed a technique 
for solving the nonlinear stream function equation for transonic flow. The 
great advantage of using the stream function, ijj, rather than the potential 
function, cj), is that rotational flows can be analyzed with tp, whereas rota- 
tion is excluded by the potential. The conservative equation satisfied by 
the stream function for rotational flow is: 



where pu = ip.,, pv = - 
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and the entropy S and the rotation w are functions of alone. Some of 
the earliest attempts at computational transonics used the stream function 
formulation, but they generally failed because of difficulties in extracting 
the density p as a function of the derivatives and The density 

is a double- valued function of ip, whereas it is single-valued in terms of 
the speed, as in equation (4). In the old formulations, the speed components 
were calculated by u = 'Py/p» v = - ip^/p, making equation (4) a nonlinear 

equation for p with a branch at the sonic speed. Hafez' idea is to obtain 
the speed components by integrating equation (3) after each iteration on 
the approach has worked nicely for irrotational flow, u e 0 everywhere, 
agreeing well with potential calculations. Work will proceed during summer 
1980 to carry out calculations for rotational transonic flow past cylinders 
and airfoils, comparing with Euler equation solutions of Salas (see this 
compendium) and potential calculations. 
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DEVELOPMENT OF A SMALL-DISTURBANCE EMBEDDED-GRID CODE 
FOR TRANSONIC WING/BODY/NACELLE/PYLON 
CONFIGURATIONS 

Perry A. Newman* 


Extension of transonic-flow computational techniques from 2-D airfoil 
sections to more realistic 3-D configurations is hindered by: (a) The compu- 

tational resources required for code development and validation as well as 
(b) the need for a suitable coordinate system. A computational grid embedding 
technique (ref. 1) solves, to some extent, both of these problems. (See 
fig. 1.) This technique was initiated by Charles W. Boppe while an NASA- 
Grumman Research Associate at Langley; continuing development has been done 
at Grumman under contract NASl-14732 (refs. 2-4). A wing-body code has been 
delivered and documented (ref. 5) while the wing/body/pod/pylon/winglet code 
(WBPPW) is in the process of being delivered to NASA. 

Comparisons of initial calculations from the WBPPW code with experimental 
results for a number of configurations are given in references 3 and 4; the 
attached figure 2 is an example. Throughout the development of this technique 
for complex configurations, simplicity has been stressed. Coordinate systems 
are essentially rectangular in character, boundary conditions are planar, 
linear interpolations are used for multiple grid interactions, a simple 
two-dimensional strip boundary-layer analysis provides viscous effects, and, 
finally a fast, easy to use fuselage modeling system yields arbitrary body 
surface normals. When compared to other computational methodology, this 
approach lacks sophistication. However, this simplicity plays an important 
role in providing the flexibility required to treat complex configurations. 
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VALIDATION OF A SMALL-DISTURBANCE 


EMBEDDED-GRID CODE FOR TRANSONIC 
WING/BODY/NACELLE/PYLON CONFIGURATIONS 
Perry A, Newman* 


Validation of the Boppe/Grumman 3-D transonic analyses codes (ref. 1) is 
being done using the experimental data base for the NASA Supercritical-Wing 
EET model. This work was initiated by E. G. Waggoner while a NASA-Vought 
Research Associate at Langley. Some results obtained from the Wing/Body code 
(WIBCO) studies are presented in reference 2. Validation studies using the 
Wing/Body/Pod/Pylon/Winglet code (WBPPW) are being continued by Waggoner under 
contract NASl-16055. 

These validations consist of comparing the results predicted by the codes 
with the carefully documented data base generated as the NASA supercritical- 
wing EET model evolved through design changes made and tested in the Langley 
8-Foot Transonic tunnel. These include addition of nacelles, pylons, and 
winglets to the clean wing-body model. The attached figure (from ref. 2) shows 
wing planform and inboard airfoil section changes made to the baseline wing- 
body configuration. The table presents a comparison of the resulting experi- 
mental and calculated force and moment increments. These comparisons show 
that the computational results are sensitive to subtle design modifications 
for the wing-body transport configurations indicating that the code can be 
used as an effective guide in the design process. Validation studies have been 
initiated on the Wing/Body/Pod/Pylon/Winglet configurations using the newly 
developed WBPPW code. 

References : 

1. See preceding compendium summary entitled "Development of a Small-Distur- 
bance Embedded-Grid Code for Transonic Wing/Body/Nacelle/Pylon Configura- 
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METHOD FOR NUMERICAL DESIGN OF A CONTOURED 


WIND-TUNNEL LINER FOR TEST OF A 

LAMINAR-FLOW-CONTROL SYSTEM FOR A 

YAWED SUPERCRITICAL AIRFOIL MODEL 

Perry A. Newman*, and 
E. Clay Anderson* 


The present numerical design procedure was developed in order to determine 
the shape of a contoured, nonporous, wind-tunnel liner for use in the Langley 
8-Foot Transonic tunnel test of a large-chord, laminar flow control (LFC), 
swept-wing panel which has a supercritical airfoil section. The procedure is 
based upon a simple idea and several existing computer codes which make it 
feasible. Basically, one determines bounding streamlines in the desired 
unbounded flow, makes all required blockage corrections, fairs into the existing 
tunnel, and finally builds the resulting contoured liner. The procedure has a 
general utility which is restricted by our current ability to calculate the 
viscous transonic flow field about arbitrary configurations. A brief descrip- 
tion and an application to the 2-D streamlined tunnel problem were given in 
reference 1. Aerodynamic lines for the design point test condition have been 
generated using the present method; data files of them are now being processed 
into tapes for numerically controlled milling of the liner blocks. The ultimate 
test of the propriety of the design procedure is the LFC experiment, currently 
scheduled to begin in mid CY 1981. 

The attached figure (taken from ref. 1) shows the agreement between cal- 
culated and experimentally determined wall shape for a 2-D streamlined tunnel 
test conducted by R. W. Barnwell and J. L. Everhart of NASA,(STAD) LRC (data 
unpublished). The symbols represent their experimentally determined wall shape; 
the curves were obtained from the present numerical procedure. 

Reference ; 

1. Newman, P. A., and Anderson, E. C. "Analytical Design of a Contoured Wind- 
Tunnel Liner for Supercritical Testing", NASA Conference Publication 2045, 
Advanced Technology Airfoil Research, Vol. 1, pp. 499-509, 1979. 
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DEVELOPMENT OF A COUPLED FINITE -VOLUME /PANEL CODE 


FOR ANALYSIS OF LFC 
SWEPT WING TESTS 
Perry A. Newman* 


Design of the wind-tunnel liner for the LFC supercritical wing experiment 
in the 8-Foot transonic tunnel (ref. 1) is for a single-test point condition. 

In order to assess and/or correct for the influence of the liner at off-design 
transonic-flow conditions, a 3-D analysis code is required. This WINGTEST 
code is being developed by Flow Research Company under contract NASl-16170. 

The test configuration to be modeled by the WINGTEST code is a constant- 
chord, variable airfoil section wing which spans a contoured-wall , transonic 
wind tunnel at an arbitrary yaw angle. This will be accomplished by modifying 
the newly developed FLO 29 full potential code to treat the test section, 
adapting a 3-D panel method code to treat the contraction section, and devel- 
oping an iterative procedure to couple the codes to obtain an interactive 
solution. This is a first-cut at arbitrary 3-D wind-tunnel simulation; if 
successful, many generalizations will open for investigation. 

Reference : 

1. See preceding compendium summary entitled "Method for Numerical Design of a 
Contoured Wind-Tunnel Liner for Test of a Laminar-Flow-Control System for 
a Yawed Supercritical Airfoil Model.'* 
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MULTIPLE LEVEL TECHNIQUES FOR COMPUTATIONAL 
FLUID DYNAMICS 

D. R. McCarthy* and R. C. Swanson** 


In recent times significant progress has been made in the application of 
the multi-grid technique in the solution of elliptic partial differential 
equations. This method attempts to eliminate relaxation stall by solving 
an equation on a sequence of course and fine grids. It also contains the 
structure necessary to relate selective local mesh refinements to the re- 
laxation process on coarser global grids. In the present work the objective 
is to demonstrate the effectiveness of the multi -grid process in transonic 
flow analysis. 

The Multiple Level Adaptive Technique (MLAT) or multi-grid procedure 
has been incorporated into the serial 3-D transonic potential flow code of 
Reyhner (Ref. 1). This modified program has been applied to the solution 
of the flow about the axisymmetric NASA TM X-2937 inlet. Through this 
application the many necessary criteria for operation of MLAT are being in- 
vestigated. The code provides four levels of mesh, the finest of which is 
161 (axial; by 81 (radial) by 5 (angular). After about one hour of CPU time, 
encompassing work equivalent to about 475 sweeps of the field, the unmodified 
program still retains residual (average change in solution at each field 
point) on the order of 10-5, and is effectively stalled. A recent MLAT run 
achieved convergence to within a residual of 10-7, with relaxation work 
equivalent to 40 sweeps. This is a substantial performance gain for the code. 
Additional cases with MLAT are now being performed. Also, the MLAT is being 
programmed into the recently completed vector 3-D transonic potential flow 
code of Reyhner. 

Reference: 

Reyhner, T. G., "Transonic Potential Flow Around Axisymmetric Inlets and Bodies 
at Angle of Attack," AIAA Paper 77-145, January 1977. 
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ALGORITHM FOR SOLVING THE FULL TRANSONIC POTENTIAL 


FLOW EQUATIONS FOR A NACELLE/WING/PYLON 
GEOMETRY USING THE CDC CYBER-203 COMPUTER 
T. A. Reyhner* and D. E. Reubush** 


During the past decade there has been tremendous progress in the field 
of transonic potential -flow computation. This work which started with small 
disturbance theory and simple two-dimensional configurations has been ex- 
tended to the full potential equation and several classes of three-dimensional 
geometries. Most extensions to three-dimensions have been based on some mapp- 
ing from the physical coordinates to a more convenient computational coor- 
dinate system. If such mappings can be found, the computational problem ex- 
clusive of the generation of the mapping, can be greatly simplified. The dif- 
ficulty with this approach is that as the geometry of objects to be considered 
becomes more complex, it becomes more difficult to find convenient mappings 
that will provide reasonable coordinates for the computation. A further 
problem is that with a complex mapping it can be very difficult to have any 
physical understanding of the solution process. 

A project has been in progress for several years to obtain solutions of 
the transonic potential flow equations in simple physical coordinates. This 
type of solution method does not require the finding of a new mapping for 
each new family of geometries. It also has the advantage that the entire 
solution procedure is in terms of physical coordinates and variables, thus any 
failures are much easier to interpret. The first work on this concept was for 
axisymmetric flow (1). Later the initial work was extended to axisymmetric 
geometries at angle of attack (2). The current work has extended the concept 
to general three-dimensional inlet, body and duct geometries and has vector- 
ized the code for use on the CDC Cyber 200 series of computers (STAR). Future 
work will incorporate wing and pylon geometries in addition to the nacelle. 

The work has been primarily done by T. A. Reyhner and associates of the 
Boeing Commercial Airplane Company with computer time and consultation pro- 
vided by NASA-Langley. 

A comparison between theory and experiment has been made for an asym- 
metric inlet tested in the NASA-Ames 40-by 80-foot tunnel. The model had a 
fan face diameter of 1.397 meters and was used with a high-bypass turbofan 
engine. Figure 1 presents a graphical display of the inlet and spinner geom- 
etry. The geometry description used in the analysis is smooth and the entire 
inlet was specified. Only half of the inlet is shown in figure 1 and curves 
are represented using straight line segments due to limitations of the plott- 
ing equipment. The comparison between theory and experiment for a run at 60° 
angle of attack and 38.9 m/s forward velocity is shown in figure 2. Agreement 
between theory and experiment is excellent. 
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Figure 1.- GRAPHICAL DISPLAY OF INLET GEOMETRY 
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A SMALL DISTURBANCE TECHNIQUE FOR PREDICTING THE 
EFFECTS OF JET EXHAUST FLOW ON THE CHARACTERISTICS OF 
TRANSPORT AIRPLANES 
Chen Sun* and William K. Abeyounis** 


One of the most complex and demanding problems during the design of tran- 
sport aircraft is the integration of the propulsion system. This problem is 
particularly difficult at transonic speeds where jet interference and entrain- 
ment effects can be significant. Current methods for predicting this complex 
flow for three-dimensional transport configurations have not been adequate. 

The purpose of this work is to extend a transonic small disturbance technique 
(Ref. 1) to more accurately account for jet interference and entrainment 
effects on transport configurations. 

The first step of this small disturbance technique is to "calibrate" the 
small disturbance code. Pressure distributions on the wing are calculated for 
a wing-body configuration using an exact transonic code. The wing-body case 
is then computed using the small disturbance code, modifying the wing geometry 
until the pressure distributions match those of the exact transonic code. The 
small -disturbance approximation is compensated for by tailoring the leading 
edge slopes and by adding correction ramps on the aft portion of the wing. 

The next step of the small disturbance technique is to make a more de- 
tailed nacelle or nacelle-pylon calculation. The present work will upgrade 
this calculation by incorporating the effects of a powered jet in a quasi- 
cylindrical fashion. The effects of jet entrainment will also be incorporated 
using a simple, isobaric mixing model. This calculation provides streamline 
slopes on a control surface surrounding the nacelle. 

The final step of the small disturbance technique is to incorporate the 
effects of the powered jet on the total configuration. This is accomplished 
by again computing the wing-body case using the small disturbance code. In 
this wing-body calculation the nacelle effects are included by imposing the 
streamline slopes from the detailed nacelle calculation on the control surface 
surrounding the nacelle. 

Reference: 

Gould, D. ; Sun, C. ; and Yoshihara, H: AFTI-III - Experiments and Calculations 

at Transonic Conditions: Symposium on Transonic Aircraft Technology, Aug. 

15-17, 1978. (Lancaster, Calif.) 
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Unsteady Flow 


Time-accurate methods designed primarily for the unsteady problem. 
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DEVELOPMENT OF FREQUENCY PLANE PERTURBATION 
METHOD FOR TRANSONIC UNSTEADY FLOWS 

Robert M. Bennett* 


One of the important and currently unsolved problems in aeronautics is 
accurate calculation of the reduction or dip in flutter speed that occurs at 
transonic Mach numbers. As one approach to addressing this problem, the 
Boeing Airplane Company has been working for several years under contract to 
the Langley Research Center to develop a transonic perturbation method for 
calculating oscillatory aerodynamics for use in conventional flutter analysis 
and design programs. The incremental forces for an oscillating wing or 
airfoil are calculated by a finite difference method as a perturbation from 
a calculated nonlinear flow field containing shock waves and other transonic 
effects. In the past, this method has failed to give meaningful results 
except at very low frequencies of oscillation. These frequency limits were 
a result of the numerical solution procedure and precluded results for normal 
flutter frequencies of oscillation. Recently the method was modified 
(under contract NASI -151 28) and good results have been obtained at 
frequencies of up to ten times the previously encountered limits and including 
the flutter range. The modifications include: (1) important revisions to 

the finite difference grid network, (2) a specially-developed, out-of-core, 
sparse matrix solver, and (3) revised treatment of boundary conditions. 

An oscillating flat plate is used as a test case for the finite 
difference method as exact solutions are available. A sample result for the 
flat plate oscillating at a reduced frequency, k = wc/2v, of 0.6 is shown 
in the figure. The real (or in-phase) and imaginary (or out-of-phase) 
components of lifting pressure are shown plotted versus the airfoil chord. 

There is good agreement between the numerical solution (symbols) and the 
exact solution (lines). Previously, this example could not be calculated 
satisfactorily as the reduced frequency of 0.6 is about ten times larger 
than the previous frequency limit. Although cases with supercritical flows 
at subsonic speeds have not yet been treated with the revised method, examples 
for a circular arc airfoil at low supersonic Mach numbers with a, detached 
bow shock wave have been satisfactorily calculated. The work is continuing 
in an effort to exploit this development by extending the method to treat 
three dimensional wings with thickness and by applying it to the calculation 
of transonic flutter boundaries for two and three dimensional wings. 


*Aeroelasticity Branch, 505-02-23, 804-827-2661 
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FINITE DIFFERENCE SOLUTION FOR 
FLAT PLATE OSCILLATING IN PITCH 

n =0.9 k = 0.6 


Frequency Plane Perturbation Method (Boeing) 
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DEVELOPMENT OF A FINITE ELEMENT METHOD 
FOR TRANSONIC UNSTEADY POTENTIAL FLOW 
Kenneth R. Kimble* 


The finite element method, developed originally for structural problems 
has many features which recommend its use for computing fluid flow for 
geometrically complex configurations. As illustrated in figure 1, the finite 
element mesh can be tailored to resolve flow singularities easily for simple 
problems. The generation of a corresponding grid for a complete aircraft 
configuration is relatively simple since the topology is almost unlimited in 
contrast to the rectangle-like mesh requirements of finite difference methods. 

In order to utilize finite elements in transonic flow an analog of 
upstream and mixed differencing had to be developed. This was done recently, 
and a comparison case is illustrated in figure 2 where the TSFOIL solution 
due to Murman is compared to our mixed flow finite element code. 

The unsteady transonic flow equation requires that the mesh be fine 
enough to resolve at least one wavelength in order to obtain accurate 
solutions. In order to perform the solution of the resulting equations (for 
which current iterative methods have been unsuccessful) an out-of-core 
equation solver (called a frontal elimination algorithm) is being adapted 
for use with the transonic finite element code. An example of an unsteady 
flow calculation for an airfoil with oscillating flap is given in figure 3. 

Current work is proceeding on the solution of the full potential 
inviscid steady and unsteady transonic equations using a finite element 
analog of the artificial compressibility method (ref. 1). 

References: 
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AN EVALUATION OF TWO-DIMENSIONAL SHOCK COMPATIBILITY 
CONDITIONS FOR LINEAR UNSTEADY TRANSONIC FLOW PERTURBATIONS 

★ 

Woodrow Whitlow, Jr. 


An important assumption in treating unsteady transonic flow problems as 
linear perturbations about a nonlinear steady state is that the aerodynamic 
loads vary linearly with the amplitude of the unsteady motion. Inherent in 
that assumption is the requirement that the shock excursion amplitude, an 
example of which is shown in Figure 1, remain sufficiently small such that 
aerodynamic loads show the necessary behavior with motion amplitude. Various 
factors, such as airfoil geometry and the amplitude and frequency of its 
motions, affect the unsteady shock movement. However, the flow conditions 
across shock waves and their movements may be determined by imposing a compati- 
bility condition at the steady state shock locations. 

The early perturbation method of Weatherill et al.^ neglected the 
unsteady shock conditions but is now being modified to include a normal shock 
condition. Cunningham^ has also derived a normal shock compatibility 
condition, and Williams^ has presented a condition that includes the effects 
of local shock inclination. Part of the current effort is also aimed at 
determining an alternate shock condition that includes the effects of local 
shock slope. 

In this study, to be conducted at Langley Research Center, the various 
compatibility conditions are to be incorporated into finite difference solu- 
tions of the perturbation equations. Where possible, comparisons with experi- 
mental data will be made. It is hoped that the results will determine the 
importance of including the effects of shock inclination in the compatibility 
conditions. Also, the aerodynamic loads will be monitored for various types 
of unsteady motions, reduced frequencies, and amplitudes to determine the 
conditions for which the assumption of linear unsteady perturbations is no 
longer valid. 
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Figure 1. Typical unsteady shock motion. 


A FINITE ELEMENT METHOD FOR THE 
PERIODICALLY OSCILLATING THIN AIRFOIL 
G. Fix, M. D. Gunzburger, R. A. Nicolaides and C. Cox* 


This investigation deals with the application of a least squares finite 
element method to the oscillating thin airfoil problem. The governing 
equation is the (small disturbance) unsteady potential transonic flow equation. 
The dependence on time is removed by assuming periodic solutions with frequency 
u), which is the frequency of oscillation of the airfoil, and the equations are 
linearized by assuming that the oscillations are small. This equation is 
recast into a first order system by introducing the velocity vector ^=(u,v) 
which is, of course, Vcf), being the potential. All boundary conditions 
are also cast in terms of the variables y^, and (}). In particular, on the 
wing we specify v and across the wake behind the trailing edge we require 
that V and the pressure, i.e. (u+iwcj)) be continuous. The domain in which 
our governing equations are valid is exterior to the airfoil. In our computa- 
tional scheme we must, of course, choose a finite computational region. On 
the outer boundary of our computational region, we replace the Sommerfeld 
radiation conditions by high order approximations to the condition [1]. 
Finally, we also insure that the Kutta condition at the trailing edge is 
satisfied. The finite element method employed consists of minimizing a suit- 
able quadratic functional J(u,v,c})) over a finite element space, which we 
choose to be a space of piecewise linear functions on selected grids. With 
proper choices of grids (see [2j ) , the method is second order accurate. Due 
to either the first order formulation or the use of a least squares principle, 
the method has the following advantages: 

1. The associated matrix problem is Hermitian and positive definite regardless 
of the size of w or whether the basic flow is or is not transonic. 

2. The achievement of second order accuracy is also independent of the above 
considerations. 

A computer code implementing the method has been developed and is 
presently in the final stages of testing. 

References: 
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IMPROVED SONIC-BOX COMPUTER PROGRAMS FOR CALCULATING 
TRANSONIC AERODYNAMIC LOADS ON OSCILLATING WINGS WITH THICKNESS 


Song Y. Ruo* 

E, Carson Yates, Jr.** 

Jerome G. Theisen* 

The potential equation for steady transonic flow is essentially nonlinear. 
In contrast, the corresponding smal 1 -perturbation equation for unsteady flow 
can be reduced to a linear equation with constant coefficients if the 
frequency of oscillation is high enough. Unfortunately, the linearization 
requires frequencies above the range of usual interest for lifting-surface 
flutter. Moreover, an important consequence of the linearization is the 
suppression of deviations are important in the propagation of disturbances 
over the lifting surface, improvement in the theory may be accomplished by 
retaining at least the effect of finite wing thickness on variations in 
mean local Mach number. This is achieved by treating steady-flow properties 
as local parameters in the unsteady flow. This assumption, equivalent to 
local linearization, permits the nonlinear differential equation for the 
velocity potential to be reduced to a linear equation with variable 
coefficients containing the local Mach number. After a nonuniform coordinate 
transformation, the equation becomes identical to the linearized transonic 
unsteady-flow equation with constant coefficients. Hence the problem can be 
solved in the transformed space by any method that is suitable for the linear 
equation, such as the sonic kernal function of the sonic-box method. 

For exploratory use in flutter analysis, the present method has been 
implemented as a modification to the sonic-box program. To illustrate its 
application, the accompanying figure shows the detrimental effect on flutter 
speed of increasing the thickness of a delta wing with elliptical cross 
section in the lateral direction. Increasing thickness from 0 to 4 percent 
reduces flutter speed by 15 percent. 

An alternate modification of the sonic-box program is also available. 

For this modification, the coordinate transformation is not used. Instead, 
the velocity potential is formulated directly in terms of local Mach number. 


*Lockheed-Georgia Company 
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EFFECT OF THICKNESS ON TRANSONIC FLUTTER SPEED 
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5. Inviscid Euler Equations 


Methods for solving the equations of inviscid, rotational flow in primitive- 
variable (p, p, u, V, w) form. Validation of codes included. 
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UPSTREAM DIFFERENCING FOR COMPRESSIBLE FLOW 


Bram van Leer* 


The one-dimensional Euler equations of compressible flow are usually 
written in conservation form when shock-wave propagation is considered and in 
diagonal or characteristic form when sound-wave propagation is considered. The 
physical content of the characteristic equations is best preserved if upstream 
differencing is used. Shock waves are properly accounted for if the scheme is 
cast in the conservation form. Van Leer has shown how to combine upstream 
differencing and conservative differencing in explicit nodal-point schemes [1] 
and control -volume schemes [2], [3]. 

In the latter scheme the characteristic equations are used to evaluate 
the fluxes of conserved quantities at cell interfaces while the conservation 
laws are used to update the cell integrals of conserved quantities. In each mesh 
a set of monotone basis functions is used to approximate the initial-value dis- 
tributions, thus eliminating numerical oscillations and providing sub-grid reso- 
lution. 

The present work aims at the following extensions of these schemes: 

a) Switching to an implicit mode wherever the local Courant number ex- 
ceeds unity (see Fig. 1). 

b) Using discontinuous basis functions for achieving the equivalent of 
shock fitting (see Fig. 2). 

c) Developing multi-dimensional forms, factorized or not, that account 
for the underlying bi -characteristic equations. 

d) Deriving efficient nodal-point versions, using the technique of flux- 
vector splitting developed by Steger, Beam and Warming [4]. 

References : 
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x/ax 1 

Figure 1 : Time-space diagram for control -volume scheme. Characteristics are 

drawn backward in timeQfrom cell boundaries. The domain of depen- 
dence of point A at t^ is included in the adjacent meshes. For 
B the domain extends beyond the left-hand mesh (Courant number >1); 
boundary data at x/Ax = -1, t > t^ may be introduced. This 
causes a dependence on C and makes the scheme implicit. 



Figure 2 : Ini tial -val ue representation. Left: Q is discontinuous in mesh 

(0,1). Right: approximation by linear basis functions. Thin 

line: least-squares approximation. Thick line: nearest monotone 

representation with same mesh integral. 

If discontinuous basis functions were included, the original 
initial-value distribution could be retrieved from mesh integral 
and mesh-boundary data. 
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FINITE DIFFERENCE SOLUTION OF A GALACTIC FLOW PROBLEM 


B. van Leer*, T. Zang, M. Y. Hussaini* 
G. D. van Albada, and W. W. Roberts 


Cosmic gas flows are distinguished from terrestrial and laboratory flows 
in deriving their main features not from shaped, solid boundaries but rather 
from gravitational forces and radiative processes. A cosmic flow problem, when 
solved numerically, may bring out lesser known properties of the finite- 
difference scheme employed. One simple, yet representative, cosmic flow problem 
has been selected for testing a number of schemes currently used in aerodynamics 
and astrophysics. 

The problem considered models, in one dimension, the periodic flow of 
isothermal galactic gas across the spiral arms of a spiral galaxy, driven by the 
gravitational field of the underlying mass distribution in the system, [1]. The 
rise and fall of the gravitational potential has the same effect as the widen- 
ing and narrowing of the channel in channel flow. Flow-dependent Coriolis forces 
add to the complexity of the flow. 

The problem has a stationary solution that may or may not include a shock, 
depending on the amplitude of the potential oscillation. The following schemes 
have been or will be tested on their ability to approximate this stationary 
solution. 


- The "beam scheme" [2] used only in astrophysics 

- Godunov's method [3] in several forms 

- Lax-Wendroff scheme with higher-order filter 

- MacCormack scheme with higher-order filter (see Fig. 1) 

- Various forms of SHASTA-based Flux Corrected Transport [4], no con- 

vergence to asymptotic solution 

- MUSCL (monotone higher-order upsteam-biased differencing)[5] 

- Pseudo-spectral method (see T. Zang and M. Y. Hussaini, "Mixed Spectral 

Finite Difference Methods for Compressible Navier-Stokes Equations," 
in this Compendium) 
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COMPUTATION OF THE SHAPE OF SLENDER STREAMS OF FLUID 
James F. Geer and John C. Strikwerda* 


The shape of a slender jet of fluid falling steadily under the force of 
gravity is studied. The problem is formulated as a nonlinear free boundary 
value problem for the potential. The use of perturbation expansions results 
in a system of equations that can be solved efficiently by integrating them 
in the stream-wise direction. Computations have been made for jets issuing 
from orifices of various shapes including an ellipse, a rectangle, and an 
equilateral triangle. Figures 1 and 2 show the cross-sectional shapes for 
vertically falling jets issuing from orifices in the shape of a triangle and 
a rectangle. In each case the jet forms sheets of fluid which propagate 
outward in directions perpendicular to the sides of the original shape. 
Surface tension effects have been ignored in these examples. 

Work has been completed for vertically falling jets without surface 
tension (Ref. 1 and 2). Current research is on nonvertical jets and jets 
with surface tension. 

The method is stable and second-order accurate (Ref. 3). This method 
appears to be new and may be applicable to other nonlinear free boundary 
value problems. 

References : 

1. Geer, J. F. and Strikwerda, J. C. "Vertical Slender Jets" Journal of 
Fluid Mechanics, to appear. 

2. Geer, J. F. and Strikwerda, J. C. "Vertical Slender Jets" ICASE Report 
79-23. 

3. Strikwerda, J. C. and Geer, J. F. "A Numerical Method for Computing the 
Shape of a Vertical Slender Jet" ICASE Report 80-7. 
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STUDY OF TRANSONIC FLOW OVER A CIRCULAR CYLINDER USING 


THE EULER EQUATIONS 
Manual D. Salas* 


During the last decade, the transonic flow over two-dimensional airfoils 
has been investigated in great detail through numerical solutions of the full- 
potential equation. It is well known that the accuracy of these solutions can 
deteriorate if the flows are strongly rotational. The purpose of this work is 
to develop benchmark solutions using the Euler equations in order to better 
assess the importance of rotational effects and the range of validity of the 
potential solutions. As a preliminary effort on the development of a code for 
general two-dimensional airfoils, a study is made of the flow over a circular 
cyl inder. 

To eliminate ambiguities in imposing the free-stream boundary conditions, 
the infinite domain outside the unit circle is mapped into the finite domain 
inside the unit circle by means of an inverse radial transformation. The time- 
dependent Euler equations are written in their A-form (ref. 1). The Euler 
equations are then integrated in time using the Gabutti finite difference 
scheme until a steady state is reached. The Gabutti scheme has the advantage 
of being stable up to Courant number 2. 

A comparison of the surface Mach number calculated using the Euler equa- 
tions and the potential equation (in nonconservative) form is shown in figure 1 
for a free-stream Mach number of 0.45. There is good agreement between the tv/o 
solutions in the forward section of the circle, with some disagreement in the 
rear of the circle. At this point in the study, it is premature to account for 
the discrepancies. The Euler solution shows a strong isentropic recompressing 
ahead of the shock which is being treated as a fitted discontinuity in the 
numerical analysis. 

The isobars for the Euler solution for the same case are shown in fig- 
ure 2. The pressure levels shown are nondimensionalized with respect to the 
free-stream pressure. 

References: 

1. Moretti, G. "The A Scheme" Computers and Fluids, Vol. 7, pp 191-205, 1979. 
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Figure 1. Calculated surface Mach number for a circle at M^ = 0.45 
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STUDY OF FLOWFIELDS ABOUT ASYMMETRIC EXTERNAL CORNERS 


Manuel D. Salas* 


A numerical study of flowfields about asymmetrical conical corners 
formed by the juncture of swept compressive wedges, which are typical of 
delta wings and rectangular inlets, see figure 1, is being conducted. 

The region of interest is the region bounded by the outer bow shock, 
the two crossflow sonic lines, and the surface of the body. This region 
is mapped into a rectangular computational plane where grid points are 
evenly distributed. The resulting mesh in the physical plane is shown on 
figure 2. 

In this region the time-dependent conical Euler equations are inte- 
grated in time using MacCormack ' s finite difference scheme, with the bow 
shock fitted as a discontinuity, until a steady state is reached. 

The main thrust of the study is an investigation of how the cross- 
flow transitions from a symmetrical configuration, where the corner is a 
cross-flow stagnation point, to an asymmetric configuration, where the 
cross flow spills over the corner as shown in figure 3. 

The details of the investigation will appear in reference 1. 
References : 

1. Salas, M. D. "A Careful Numerical Study of Flowfields About Asymmetric 
External Conical Corners" AIAA Paper No. 80-1329, AIAA 13th Fluid and 
Plasma Dynamics Conference, July 1980. 
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INTERACTION OF A TWO-DIMENSIONAL SHOCK WAVE WITH A VORTEX 
Manuel D. Salas and Sui-Kwong P. Pao 


This investigation considers the time-dependent interaction of a vortex 
field and a shock wave inside a shock tube. The numerical investigation 
attempts to duplicate the experimental study reported in reference 1. The 
primary interest in this phenomenon concerns the noise generated by the shock- 
vortex interaction. The phenomenon can also be important in learning how to 
handle numerically the shock-vortex interaction that occurs in high-speed 
fighter wings. 

In the numerical study, the two-dimensional Euler equations are integrated 
in time using the MacCormack finite difference scheme to simulate the transient 
motion of a shock wave over a given vortex field. The region of interest, 
see figure 1, is the region bounded by the shock-tube walls, the shock wave 
on the right, and the given supersonic inflow conditions at x = 0 on the 
left. 


As the initially straight shock moves through the vortex field, the shock 
profile deforms into an S shape. As the interaction continues a kink develops 
in the shock profile, where a secondary shock wave forms. See figure 2. 

These numerical results are in agreement with the experimental observations 
of reference 1 . 

The appearance of the secondary shock wave complicates the numerical 
calculation of the shock making the MacCormack scheme inappropriate for this 
problem. We presently plan to replace the MacCormack scheme with the X-scheme, 
see reference 2, to alleviate this problem. 

References: 

1. Naumann, A.; and Hermanns, E. "On the Interaction Between a Shock Wave and 
a Vortex Field" AGARD-CP-131 . 

2. Moretti , G. "The X Scheme" Computers and Fluids, Vol . 7, pp. 191-205, 1979. 
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Figure 1.- Sketch of configuration being considered. 



Figure 2.- Isobars for shock-vortex interaction. 



CALCULATION OF SUPERSONIC FLOW PAST AIRCRAFT WITH SHOCKS 
Manuel D. Salas* and James C. Townsend** 


Under a grant to the Polytechnic Institute of New York, a computer 
code is being developed to solve the three-dimensional, supersonic, steady 
Euler equations for the flow field about an arrow wing aircraft. The hyper- 
bolic character of the equations allows integration of the equations from 
one axial station to the next. At each axial station a series of conformal 
mappings is used to transform the region between the bow shock and the 
fuselage of the aircraft into a rectangular computational region. The 
Euler equations are then integrated using the A-scheme.^ 

A comparison between experimental and calculated surface pressures 
is shown in figure 1 at two stations. These results are at M = 2.36 and 
zero angle of attack. 

Presently, the code is being extended to include the angle-of-attack 
capabi 1 ity. 

References ; 

1. Moretti, G. "The X Scheme" Computers and Fluids, Vol. 7, pp. 191-205, 
1979. 
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DEVELOPMENT OF A MULTI-GRID METHOD 
FOR THE EULER EQUATIONS 

Stephen F. Wornom 


Time accurate and relaxation procedures for the compressible inviscid 
fluid conservation-law-equations (Euler equations) are frequently less 
efficient than procedures developed for potential flow. This is so because 
with the Euler equations one has to deal with a system of nonlinear partial 
differential equations as opposed to a single nonlinear partial differential 
equation in potential flow. The consequence of a system of equations is that 
a spectrum of associated eigenvalues must be considered so the equations are 
usually stiffer. Of more import, with a system of equations, the spectrum of 
associated eigenvalues can appear (as in subsonic flow) with both positive and 
negative real parts. This last restriction severely limits the choice of 
allowable difference operators and iterative or time-dependent procedures. 

The purpose of this study is to develop and thoroughly test a point 
relaxation algorithm for the Euler equations which is made possible by use 
of flux vector similarity splittings. In this approach, the flux vectors 
are split such that their associated Jacobian matrices have either all posi- 
tive or all negative real parts. The split vectors can then be specially 
differenced using dissipative forward or backward operators. Simple point 
iterative methods can thus be used as a solution procedure. A major part of 
the study will be devoted to improving the iterative convergence rate using 
the multi-grid relaxation approach. 

This research is being conducted by Dr. J. L. Steger of Flow Simulations, 
Inc., under contract NASl-14517. 
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SOLUTIONS TO THE EULER EQUATIONS USING A STREAMLINE- 
CHARACTERISTICS COORDINATE SYSTEM 
Lawrence Sirovich* and James C. Townsend** 


Although current finite-difference methods offer a reasonable approach 
to computational fluid mechanics for supersonic flows and can produce good 
results, they have a variety of difficulties which make it worthwhile to 
investigate new approaches. The present work, which is being done at Brown 
University under grant NS6 1617, uses flow field physics, analytical methods 
and known approximations in developing a rapid and accurate numerical code. 

The method has been applied to unsteady one-dimensional flows and to steady 
two-dimensional flows with attached shockwaves. 

Figure 1 shows such a flow over a two-dimensional airfoil. The two- 
dimensional, inviscid steady flow equations are written in characteristic 
form with entropy, flow angle and Mach angle as independent variables in 
terms of coordinates along the streamlines and principle (C'*’) characteristics. 
An approximate solution, found by assuming that the flow angle 0 is constant 
along the C+ characteristics, is used to map the entire upper half plane in- 
to a finite region in the computational plane (fig. 2). The body surface and 
the streamlines become straight horizontal lines in the mapped plane, and the 
C* characteristics become vertical lines. The curved lines in figure 2 are 
the leading-edge and trail ing-edge shock waves, and the point where they meet 
(a = 1 ) corresponds to conditions at infinity, where these shocks have weakened 
to Mach 1 i nes . 

With the leading-edge shock location fixed in the computational plane, a 
simple iteration is applied to satisfy the relations along the C~ characteris- 
tics (neglected in the approximate solution). Figures 3 and 4 compare the 
approximate solution for a typical case with the results after just three 
iterations, where the solution is converged to within 1 percent throughout 
the flow field. The converged solution (which is the solution of the exact 
equations) includes the transformation locating the shocks, streamlines and 
characteristics in the physical (x,y) plane. The usual physical flow varia- 
bles (pressure, velocity, etc.) are recovered from the computational 
variatiles through simple relations. 

Work is continuing with the eventual goal of producing a code for fully 
three-dimensional supersonic flows. 

* Brown University, 401 '-86^-21 T4 
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Figure 1.- Supersonic flow past a two- 
dimensional airfoil, showing shock 
waves (heavy lines), streamlines and C+ 
characteristics (solid lines), and C~ 
characteristics (dashed lines). 


leading-edge trail ing-edge 



Figure 2.- Approximate solution in 
computational plane for Mach number 2.5. 
Streamlines are constant a ; C'*' charac- 
teristics are constant 3. 




a 


Figure 3.- Approximate and converged Figure 4.- Approximate and converged 
solutions for flow angle 0 on each C+ solutions for Riemann invariant on each 
characteristic of figure 2. streamline of figure 2. 
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CALCULATION OF TWO-DIMENSIONAL INLET 


FLOW FIELDS IN A SUPERSONIC FREE STREAM 
Wallace C. Sawyer* 


It is of great interest to the designers of supersonic airbreathing 
missile configurations to understand the complexities of the flow field about 
the configuration forebodies and the inlet. Analysis methods are of consider- 
able importance to produce efficient designs, e.g., to maximize total pressure 
recovery, and to minimize cowl drag due to spillage at off-design conditions. 

Nielsen Engineering and Research, Inc., under contract to NASA/LaRC 
(ref. 1-4) has developed a method to calculate two-dimensional inlet flow 
fields in a supersonic free stream by an implicit, shock-capturing, finite- 
difference method. The Euler equations are subjected to a general curvilinear 
transformation and a body-fitted coordinate system is employed. The method 
is used to solve super-critical, critical, and subcritical flow fields which 
are simulated by prescribing appropriate conditions at the inlet outflow 
boundary. 

Typical results are shown in figure 1 for an external compression inlet 
composed of a 10° ramp and a 20° cowl. These Mach 2 results are in excellent 
agreement with both experimental data and results of an exact solution. 

Extensions of the procedure to the case of a nonuniform-free stream is 
straight forward and is being pursued. Work is also in progress for the 
inclusion of laminar and turbulent viscous effects. A two-dimensional inlet 
model is under construction, the test results will be compared with theory 
for validation of the computational techniques. 

References : 

1. Biringen, S. H.; and McMillan, 0. J.: Calculation of Two-Dimensional Inlet 

Flow Fields in a Supersonic Free Stream - Program Documentation and Test 
Cases. NASA CR-3221 , March 1980. 

2. Biringen, S. H.; and McMillan, 0. J.: Calculation of Two-Dimensional Inlet 

Flow Fields in a Supersonic Free Stream by an Implicit Marching Code With 
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3. Chaussee, D. S.; and McMillan, 0. J.: A Supersonic, Three-Dimensional Code 
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4. Chaussee, D. S.; and McMillan, 0. J.: A Supersonic Three-Dimensional Code 

for Flow Over Blunt Bodies - Program Documentation and Test Cases. 

NASA CR-3224, February 1980. 
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X-COORMt<ATE OF EXTERNAL SURFACES 



TNO-DIMENSIONAL SINGLE WEDGE INLET SUPERSONIC FLOW M„= 2.0 


THE FLOW RESULTING FROM THE INTERSECTION OF AN OBLIQUE SHOCK 
WAVE WITH A CYLINDER ALIGNED WITH A STEADY SUPERSONIC FLOW 

James C. Townsend* 


This study addresses the general store-carriage and separation problem. 
Current analysis methods use linearized supersonic flow theory which does not 
account for the complex shock wave interactions occurring when a store (such 
as a missile) is carried in close proximity to an aircraft. A program has 
been started to examine the details of the flow for this situation, both 
numerically and experimentally, by first breaking the general problem down 
into several simple unit type problems. 

This particular unit problem is the first reflection of a wing shock 
from a store. The wing shock is represented as a planar oblique shock wave, 
and the store it intersects is represented as an infinite cylinder aligned 
with the supersonic free-stream flow. The free-stream Mach number and the 
incident shock angle determine the flow conditions between the incident and 
reflected shock waves, so only the reflected shock conditions and the flow 
behind it need to be computed. As shown in figure 1, the incident shock is 
simply reflected where it first intersects the cylinder, but at some point as 
it wraps around the cylinder the simple reflection can no longer be sustained 
and the reflection changes to the more complex lambda configuration. 

For computation the region in each cross-section plane between the 
reflected shock and the cylinder is mapped into a rectangular region (fig. 2). 
The steady-state Euler equations are integrated using a finite-difference 
method to march in the free-stream direction, that is, along the cylinder 
axis. The problem is in the initial formulation stage with no computed 
results available yet. 

A complementary experimental study is also planned in support of this 

work. 
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COMPARISON OF RESULTS FROM A FINITE DIFFERENCE CODE 
WITH RESULTS FROM EXPERIMENTS AND OTHER METHODS 
Emma Jean Landrum* 


A number of computer codes have been developed or are being developed 
which show promise of significantly improving the estimation of aerodynamic 
characteristics for arbitrarily-shaped bodies at supersonic speeds. 

One of these codes is a shock-fitting finite difference method (ref. 1) 
based on general numerical scheme for solving the Euler equations for 
supersonic flow about complex configurations. A finite difference marching 
scheme of second order accuracy is used to compute the complete flow field 
between the body and the bow shock. A conformal mapping of the region between 
the bow shock and the body into a rectangular region produces the computational 
grid. When embedded shock waves occur, the mesh is adjusted so that the mesh 
lines coincide with the shocks and the Rankine-Hugonoit relations are satisfied 
explicitly across each of the shock waves; this is the floating, shock-fitting 
method. 

Results for an analytically-developed forebody (which resembles the cock- 
pit region of a fighter aircraft) are compared with experiment and four 
linearized theory methods in figure 1 for Mach numbers of 1.7 and 4.5 at an 
angle of attack of 5°. The finite difference method provides excellent agree- 
ment with experiment at all Mach numbers. In addition, computer solution time 
for the shock-fitting-finite-difference method is competitive with other codes 
(ref. 2). 


References: 
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(a) e = 0°. 

Fig. 1.- Analytically developed forebody, 
a = 5°. 
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Fig. 1.- Concluded. 



APPLICATION OF A SOLUTION ON THE EULER EQUATIONS TO 
THE PREDICTION OF STATIC PERFORMANCE OF 2-D C-D NOZZLES 

Mary L. Mason* 


Advanced computational models are currently being developed for the 
solution of the two-dimensional and three-dimensional Euler equations 
applied to internal nozzle flow and exhaust jet flow. These algorithms 
must be evaluated by comparisons of the theoretical results with experimental 
data. Comparisons of theory and data are essential in determining the 
accuracy of the theoretical results and in defining a valid range of appli- 
cation of the theory. 

An experiment has been conducted in the 16-foot transonic tunnel static 
test facility to define internal performance characteristics for five two- 
dimensional C-D nozzles. Internal static-pressure data (p/ptj) along 
the flaps and sidewalls and internal thrust ratio data (F/Fi) were taken 
at nozzle pressure ratios ranging from 2.0 to approximatley 9.0. The in- 
ternal performance data may be applied to the evaluation of computational 
models for prediction of nozzle internal flow behavior. 

Since the internal flow data exhibit predominantly two-dimensional char- 
acteristics, the force and pressure measurements were compared with the 
results of the two-dimensional inviscid model of Cline (ref. 1). This algo- 
rithm solves the two-dimensional Euler equations for nozzle internal flow 
and exhaust jet flow, with shock effects incorporated using an artificial 
viscosity. The results of the data-theory comparison are given in fig. 1 
for one of the experimental 2-D C-D nozzle configurations. The x coordinate 
used in the figure represents the distance from the nozzle throat, positive 
downstream. The I coordinate used to nondimensionalize x is the distance 
from the nozzle throat to the nozzle exit. Agreement of data with theoret- 
ical results is good for both thrust ratio and static-pressure data. The 
theory is limited in application to near-design conditions or underexpanded 
conditions. Since the algorithm is inviscid, it is inadequate for modeling 
severely overexpanded nozzle flow where internal separation may occur. 
However, at design nozzle pressure ratios and above, the good agreement be- 
tween data and theory indicates the effective application of two-dimensional 
inviscid theory to preliminary design-performance prediction for 2-D C-D 
nozzles. 

Reference: 

Cline, Michael C.: NAP: A Computer Program for the Computation of Two- 

Dimensional, Time-Dependenc Inviscid Nozzle Flow. LA-5984 (Contract W-74 
05-EN6-36), Los Alamos Scientific Lab., Jan. 1977. 
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A SOLUTION OF THE EULER EQUATIONS 
FOR THE FLOW IN THE ENGINE NACELLE-PYLON-WING 
REGION OF TRANSPORT AIRPLANES 
William B. Compton III* 


The integration of a propulsion system with the airframe of a transport 
airplane must be given close attention if an aerodynamical ly efficient 
configuration is desired. The problem becomes even more critical at tran- 
sonic speeds where strong compressibility effects occur and shock waves are 
possible. Currently, calculations of the flow over the relatively complex 
engine nacelle-wing-pylon geometry can only be made with theories on the 
order of fully linear potential flow coupled with a compressibility 
correction. Calculations using the Navier-Stokes equations are limited to 
the flow over less complex configurations such as the afterbody of a nacelle. 

This project was initiated to develop a more exact numerical solution of 
the flow in the region of an airplane's engine nacelle and wing. The Euler 
equations will be solved using an implicit, time-dependent, finite-difference 
approach. The Beam and Warming factored ADI scheme was chosen as the 
numerical algorithm. The solution will be good for highly compressible flows, 
and will allow for streams having different entropy levels. 

Figure 1 shows the computational domain. Within this region, a boundary- 
fitted coordinate system is generated by solving Laplace's equation, the 
Thompson-Thames-Mastin method of generating grids, at cross-sectional planes. 
The curvilinear coordinate system produced in the physical space is then 
transformed, into a cylindrical coordinate system for the computations. 

Figure 2 shows a cross-section of the coordinate grid in physical space. 

As a first step in getting the program operational, the grid shown in 
figure 2 is being used to compute the flow over a cylinder. If this compu- 
tation establishes the suitability of the coordinate system, the program will 
be extended to 3 dimensions with an engine nacelle-pylon-wing geometry. 
Initially, the airplane geometry will be simplified and only configurations 
with pointed leading edges and a vertical plane of symmetry will be considered. 
Future extensions of the program are planned to include the capability of 
handling blunt leading edges and swept wings. 
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PHYSICAL COMPUTATIONAL 

PLANE PLANE 


Figure 1. - Computational domain. 



Figure 2. - Cross-section of the coordinate 
grid in physical space. 
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6. Navier-Stokes Equations 


Methods for solving the Navier-Stokes equations. 
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Incompressible Flows 


Methods designed primarily for solution of incompressible flow. 



BOUNDARY-LAYER RESOLVING DIFFERENCE EQUATIONS FOR 
STEADY, INCOMPRESSIBLE NAVIER-STOKES EQUATIONS 
Thomas B. Gatski*, Chester Grosch** and Milton E. Rose** 


A particularly natural (spline-type) approximation basis for solving a 
boundary value problem for a linear ordinary differential equation is a patch- 
basis (ref. 1); under certain circumstances there results a set of difference 
equations whose solution exactly solves the differential equation at mesh 
points. This idea has been independently developed and applied to singular 
perturbation problems in order to demonstrate that accurate difference schemes 
can be developed which help avoid the need for mesh refinement in boundary 
1 ayers . 

This research investigates an extention of this method to the steady- 
state vorti city-stream function formulation of the Navier-Stokes equations 
in two dimensions. A central part of the study aims at comparing classes of 
computational boundary conditions which have been proposed for these equations. 
A Newton iteration scheme using a multi grid technique is being explored in 
order to enable experiments to be conducted economically. 

References : 
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SPLIT-VELOCITY NAVIER-STOKES TECHNIQUES FOR HIGH 
REYNOLDS NUMBER FLOWS 
Douglas L. Dv/oyer* 


The spl it- velocity Navier-Stokes solution procedure was originally 
developed by Dodge (ref. 1) for the solution of internal flow problems. 

The method consists of splitting the velocity vector into its rotational 
and irrotational parts in such a way that the irrotational part of the 
velocity vector is identified with pressure. This split allows one to 
solve a modified form of the potential equation for the pressure and trans- 
port equations for the rotational part of the velocity. The procedure should 
be well -suited to high Reynolds number flows since the equations for the 
rotational part of the velocity need only be solved in the viscous dominated 
parts of the flowfield. 

In the present work, the extension of the split-velocity method to 
external flows is investigated. Incompressible laminar flow over a flat 
plate (semi-infinite and finite) was chosen as a model problem. In the 
laminar case the calculated results agreed very well with asymptotic theory 
solutions for both the semi-infinite and finite cases at high Reynolds 
number. Further, computation time was found to be virtually independent 
of Reynolds number, a typical case requiring less than 1 minute of CYBER-175 
time on a 40 X 60 grid. 

A comparison of the calculated skin friction along a turbulent semi- 
infinite flat plate with a curve fit of experimental data is shown in the 
figure at a Reynolds number per foot of 10®. There is good agreement be- 
tween the present results and the data with the calculation requiring less 
than 2 minutes of CYBER-175 time, v/hich is about 30% more computer time than 
a laminar case under the same conditions and on the same grid. 

The split velocity methods appear to be quite accurate and efficient 
for calculating external high Reynolds number flows. Based on these results, 
the method is being applied to more complex flows and a compressible extension 
is being developed. 

References : 
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0 METHOD VERIFIED FOR HIGH RN 
TURBULENT FLOW. 


0 SAME GLOBAL CONVERGENCE 
HISTORY AS LAMINAR CASE. 


0 30% MORE COMPUTER TIME THAN 
LAMINAR CASE. ADDITIONAL 
TIME DUE TO EVALUATION OF 
TURBULENCE MODEL. 


0 PLANS 

• Turbulent Incompressible 
Airfoil Calculation 

• Vector Processor Algorithms 
e Compressible Extension 


INCOMPRESSI BLE HIGH RN FLAT PLATE FLOW 
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Figure 1.- Velocity split procedure for solving Navier-Stokes 
equations. 



FULLY VISCOUS MULTI-ELEMENT AIRFOIL SOLUTIONS 

★ 

Frank C. Thames 


In order to obtain more accurate solutions to multi -element airfoil 
problems, an approach is being developed which is based on the full Navier- 
Stokes equations. Key features of the method include: 

0 Use of boundary fitted coordinates (ref. 1) as illustrated in 
figure 1. 

0 Simple, algebraic turbulence models (ref. 2). 

0 Velocity-split fluid equations (ref. 3). 

0 Multi-grid solution algorithm (ref. 4). 

Boundary-fitted coordinates are essential in viscous flow solutions to ration- 
alize the nearly singular behavior of the Navier-Stokes equations in the 
neighborhood of solid boundaries. For the present simple, algebraic turbu- 
lence models will be used. More complex ones may be added later. The velocity- 
split Navier-Stokes equation formulation allows the flow problem to be split 
into its constituent potential and viscous parts and each part can be solved 
on separate grids whose spacing is appropriate to the length scale of the 
particular type of flow. The choice of the multi -grid algorithm was based 
primarily on the method's outstanding convergence rate. 

The work to date has been focused on the grid generation task and the 
evaluation of relaxation schemes to drive the multi -grid procedure. A sample 
convergence history is given in figure 2. Current results indicate that the 
multi -directional schemes work best. 
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Figure 1.- Mul ti -element airfoil computational grid. 
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factorization relaxation driver. 
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Compressible Flows 


Methods designed primarily for solution of compressible flow.^ 
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NAVIER STOKES SOLUTION FOR THE DYNAMIC STALL 

OF HELICOPTER AIRFOILS 

★ 

Warren H. Young, Jr. 


The limitations on helicopter flight envelopes caused by vibratory loads 
must be predicted in order to assess helicopter agility, blade and control 
system stresses, and ride quality. The objective of this analytical work is 
the prediction of unsteady airloads from dynamic stall. The airloads are 
necessary to predict blade vibratory response in stall flutter. The analysis 
will be used to analytically evaluate new airfoil design on the basis of 
susceptibility to dynamic stall onset and point of recovery from stall. 

The analysis is a solution of the Reynolds-averaged Navier Stokes equa- 
tions for two-dimensional, time dependent, subsonic, compressible flow 
about an airfoil. The equations are formulated in non-conservative form with 
the energy and mass fluxes as the independent variables. Total temperature 
was assumed to be constant. The turbulent shear stress has been modeled by 
an algebraic length scale in the turbulent kinetic energy equation. Flow 
fields containing laminar, transitional, and turbulent flow are calculated 
without interruption of the solution scheme. 

The finite-difference equations are written in a non-orthogonal , highly 
stretched coordinate system. The grid allows computation for airfoils of 
general shape, efficiency in minimizing the number of grid points and consid- 
erable control in the placement of grid lines. The finite-difference equa- 
tions are solved by a linearized block implicit scheme based on the Douglas- 
Gunn ADI procedure (ref. 1). The calculation scheme allows large time steps 
and is also suitable for steady solutions (ref. 2). 

This work is continuing at Scientific Research Associates, Inc., under 
contract NASl-15214. Recent results (fig. 1) for a highly stalled airfoil 
at one million Reynolds number show the important flow features that were 
measured (fig. 2) by laser velocimetry. Based on these results, modifications 
will be made to better resolve the turbulent shear in the free shear layer 
above the airfoil. Possible future work may include weak shocks in the 
solution. 

References : 
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Viscous Analysis of Tip Vortex Formation 
Warren H. Young, Jr.* 


Recent flight tests of the UH-1 helicopter have demonstrated benefits in 
noise reduction and performance optimization from modifications of tip vortex 
size. The purposes of the present program are (1) to identify the relation 
between blade tip planform and the viscous rollup of the tip vortex, (2) to 
establish benchmark calculations for simpler, more approximate models of the 
tip vortex formation, and (3) to compare the aerodynamic performance of 
proposed tip shapes. 

A prototype turbulent solution has been completed under contract NASl- 
14904. The initiation, formation, and viscous diffusion of a vortex over the 
upper surface of a thin tip has been calculated. The theory accounts for the 
vorticity in the lower surface boundary layer (Fig. 1) being converted around 
the tip, separating from the surface, and forming the tip vortex. The 
Reynolds-averaged Navier-Stokes equations, with surface boundary conditions 
and vortex stretching approximations, are solved in calculation planes per- 
pendicular to the freestream direction. A percursor inviscid solution is 
required to generate the streamwise pressure gradient as input to the present 
analysis. The present scheme is parabolized in the freestream direction to 
allow marching through successive calculation planes. The elliptic nature of 
the Navier-Stokes equations is retained in the two directions that lie in the 
calculation plane. Solution in the plane is by an ADI scheme (Ref. 1). 

Results of the prototype solution are shown in figure 2. The arrows are 
velocity vectors in a plane normal to the freestream direction. The tip 
vortex rollup, which was calculated to begin near 30% chord, has moved above 
and slightly inboard of the tip. The calculation has produced a realistic 
qualitative results and, upon refinement, will be suitable for a parametric 
comparison of helicopter rotor blade tip shapes. These results have been 
published in reference 1. 

The development of the technique is continuing under contract to 
Scientific Research Associates. The current tasks in improving the solution 
include exact satisfaction of the no-slip condition, restoration of the vortex 
stretching terms, and extension to planforms other than rectangular. Possible 
future tasks include development of a coordinate system to allow the thickness 
and curvature about the tip to influence the vortex formation. The results 
will be used to both verify and improve the more approximate, but simoler, 
percursor inviscid solution. 
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1. Shamroth, S. J., and Briley, W. R. : A Viscous Flow Analysis for the Tip 

Vortex Generation Process. NASA CR-3184, Oct. 1979. 
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Figure 2.- Velocity vectors in the calculation plane at 75% chord. 



COMPUTATION OF THE PERPENDICULAR INJECTOR 
FLOW FIELD IN A HYDROGEN FUELED SCRAMJET 
J. Philip Drummond and Elizabeth H. Weidner* 


Work is underway at the NASA Langley Research Center to develop a 
hydrogen fueled supersonic combustion ramjet (scramjet) capable of operating 
at hypersonic speeds , in the atmosphere. Hydrogen fuel is introduced into the 
engine in both a parallel and perpendicular direction to the engine primary 
flow. This approach tailors the injection to an optimum heat release schedule 
in the engine over the flight Mach number range of interest. Efficient engine 
design requires a detailed understanding of the flow field near fuel injectors. 

A computer program has been developed at Langley to begin studying the 
flow field near one or more perpendicular fuel injectors that are injecting 
hydrogen fuel at near sonic conditions into a ducted supersonic crossflow. 

The flow field is highly turbulent and contains regions of subsonic recircu- 
lating flow on either side of the injector (fig. 1). Therefore, solutions to 
the two-dimensional elliptic Navier-Stokes and one or more species equations 
are required to describe the flow. These governing equations are integrated 
by time relaxation using the time-split finite difference technique of 
MacCormack (ref. 1). 

Flow field results for a 6 cm high by 10 cm long duct into which cold H„ 
(242 K) is being injected are given in figures 2-6 (ref. 2). Note the strong 
bow shock ahead of the injector, the expansion of the primary flow over the 
injector, and the recirculating flow upstream and downstream of the injector. 
Significant penetration of the hydrogen jet into the primary flow, and con- 
vection/diffusion of H„ upstream of the jet is also apparent in figure 6. 

These results agree qualitatively with experiment. A quantitative comparison 
of program results with a data case is given in reference 2. 

Current work with the program includes addition of a chemistry model to 
study H 2 -air reaction and improvement of the turbulence model. The analysis 
of a two-dimensional inlet/ combustor with H 2 injection will then be attempted. 
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Fig. 1 Model problem for flow between two struts. 
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Fig. 2 Velocity vector field in lower half duct. 
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Fig. 3 Magnified velocity vector field about lower 
injector. 
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Fig. 4 Static pressure contours in lower half duct 
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Fig. 5 Static temperature contours in lower half 
duct. 
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Fig. 6 Hydrogen mass fraction contours in lower 
half duct. 



NUMERICAL ANALYSIS OF THE INLET FLOW OF A SCRAMJET ENGINE 

A jay Kumar* 


NASA Langley Research Center is currently engaged in developing a hydrogen- 
fueled supersonic combustion ramjet (scramjet) engine for hypersonic speeds. 

It has a fixed-geometry inlet with wedge-shaped sidewalls. Sweep of these 
sidewalls, in combination with a recess in the cowl, allows spillage to occur 
efficiently with fixed geometry of the inlet. Inlet compression is completed 
by three wedge-shaped struts located at the minimum-area section. 

The objective of this work is to develop a numerical code to analyze the 
scramjet inlet flow. In an attempt to understand some of the aspects of this 
flow, a two-dimensional model problem (fig. 1(a)) is formulated. In this pro- 
blem, the top surface of the inlet produces a six-degree compression. The cowl 
plate is located in such a way that there is some space left between the cowl 
plate and the shock from the top surface to allow for flow spillage. A four- 
degree half-angle strut is placed in the flow as shown in the figure. 

Two-dimensional N.S. equations are used in the conservative law form to 
describe the inlet flow field. Boundary-fitted curvilinear coordinates are 
used to transform the physical domain into a rectangular domain (fig. 1(b)) with 
uniform mesh spacing. The strut EFGH is transformed into a slit E'G' which is 
coincendent with one of the mesh lines. The transformed equations are solved 
by an explicit time-asymptotic, finite-difference method due to MacCormack. 

Results presented here are for the laminar flow of air with frees tream 
conditions taken as: = 5, = 1 atm, and T^ = 293 K. Figures 2 and 3 

show the plots in the region downstream of line°°zz' in figure 1(a). Upstream 
of this line, the shock from the top surface turns the flow downward by about 
six degrees. Due to this, some of the flow spills out of the inlet. Figure 2 
shows the surface pressure distributions. It is seen that the pressure on the 
top surface remains constant downstream of the compression corner until the 
shock from the leading edge of the strut's upper surface hits the top surface 
and gets reflected. This produces a large pressure increase on the top surface. 
Downstream of this point, the expansion waves from the strut's shoulder interact 
with the reflected shock and decrease the pressure on the top surface. The 
pressure on the strut's upper surface shows an increase due to the leading edge 
shock and than a decrease due to the expansion waves at the shoulder. The 
pressure on the lower surface of the strut shows two expansion waves produced 
by the strut's leading edge and the shoulder. These expansion waves interact 
with the shock from the leading edge of the cowl. 

Figure 3 shows the plot of the velocity vectors. It is seen that the 
shock from the strut separates the flow on the top surface of the inlet. A 
small separated region is also produced near the trailing edge of the strut 
which is probably caused by the cowl shock. The flow downstream of the strut 
becomes very complex due to the multiple interactions of the shocks and the 


*01d Dominion University, Norfolk, VA, under Contract NASl-15930, 804-827-2803 


142 



expansion waves as seen in figure 4 which shows pressure variation across the 
inlet at X = .1089 m, just upstream of the inlet outflow boundary. 

A mesh size of 51 x 51 is used in the above analysis. The converged 
solution is obtained in 2150 time-steps which required two minutes of computing 
time on CDC-CYBER-203 computer. 
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Fig. 1 Configuration of model 
inlet problem 
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Fig. 2 Surface pressure distributions 


Fig. 3 Velocity vector field down- 
stream of line zz' 



Fig. 4 Pressure variation across the 
inlet just upstream of outflow 
boundary 
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GIM/STAR: EXTERNAL AND INTERNAL FLOWS FOR HYPERSONIC VEHICLES 

L. W. Spradley, J. J. Stalnaker, A. W. Ratliff,* 

J. L. Hunt and J. P. Drummond** 

The General Interpolants Method (GIM) is a hybrid procedure for generating 
numerical analogs of the three-dimensional Navier-Stokes equations. A code 
using GIM is now available on the Langley CDC/STAR computer (ref. 1). The 
GIM/STAR effort is an interrelated program, being pursued by both Langley's 
Hypersonic Aerodynamics Branch (HAB) and Hypersonic Propulsion Branch (HPB) 
through coordinated contracts with Lockheed/Huntsville. 

The GIM/STAR code combines the finite element geometry point of departure 
with finite difference computational analysis. Complex 3-D flow domains can be 
modeled due to the allowance of regional numerical transformations. Grid 
planes and metric coefficients are generated in the transformed coordinates; 
computations are made in the physical plane after the application of inverse 
transformations. Steady-state solution of the governing equations is obtained 
with a transient procedure using the MacCormack unsplit finite difference 
numerical integration technique. 

Two complex flow types, relevant to Langley's hypersonic program, that have 
been computed with GIM/STAR are shown in the figures. The first illustration 
is that of a Mach 5 flow into a two-dimensional inlet, configured to tax the 
code's capabilities. The external flow was treated inviscidly. The shock 
generated by the 25° ramp angle is shown as captured by the code; 34 percent of 
the flow is spilled at the cowl lip as predicted by theory. Inside the duct, 
the Navier-Stokes equations were solved with a no-slip condition on the upper 
wall and an assumed laminar boundary layer. The turning of the internal flow 
was accomplished by a strong expansion in the upper initial portion of the 
duct diffuser and a strong shock emanating from the cowl lip. The shock im- 
pinged the internal upper wall and resulted in separation of the boundary 
layer and formulation of a recirculation zone. The second illustration shows 
the pressure contours calculated in and about a three-dimensional supersonic 
exhaust plume emersed in a hypersonic external flow. The strong pressure 
gradient outboard the nozzle lip indicates the presence of a shock created by 
the underexpanded plume spilling over into the external expanding flow. 

Work is currently underway to enhance the overall computational speed of 
the code in supersonic viscous and inviscid flow domains by incorporating a 
forward-marching scheme in combination with transient relaxation at each cross 
plane. This algorithm allows direct solution of a parabolic form of the 
governing equations without the necessity of a pressure iteration and is 
readily adaptable to a hyperbolic solution of the Euler equations as well. 
Turbulence and chemistry models are to be incorporated into the code during 
the next 12 months. 
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MIXED SPECTRAL/FINITE DIFFERENCE 
METHOD FOR COMPRESSIBLE NAVIER-STOKES EQUATIONS 
Thomas A. Zang and M. Y. Hussaini* 


In high Reynolds number flows the numerical dissipation resulting from the 
truncation error of the advection terms may overwhelm the physical dissipation 
from the viscous terms. A numerical scheme for compressible, mildly viscous 
flows which employs a standard finite difference approximation for the viscous 
tei^ but resorts to a pseudo-spectral method for the advection terms is ex- 
amined here. This method offers the prospect of greatly reducing the major 
source of numerical error without dramatically increasing the computational 
effort (unlike a fully pseudo-spectral method which is quite costly to apply to 
the complicated compressible viscous terms). 

Several model problems have been considered. The first is the familiar 
Burger's equation with periodic boundary conditions. A spectral/finite-differ- 
ence (SFD) method was compared with a second-order, fully finite difference 
(FFD) method. Explicit time-differencing was applied to the advection term 
and implicit time-differencing to the diffusion term. Figure 1 illustrates 
the results. The solid curve is the exact solution and the open circles are 
the numerical results. Not only does the SFD method produce high accuracy away 
from the shock, but it produces a fairly narrow shock as well. 

The SFD method has also been applied to a set of viscous, compressible, 
quasi-one-dimensional euqations that approximate the flow in a nozzle of slowly 
varying cross section. Here, the pseudo-spectral portion uses Chebyshev poly- 
nomials. The steady-state numerical solution for one case is shown in Figure 
2. Two Chebyshev grids have been patched together at the throat in order to 
improve the resolution at the sonic point. Some smoothing has been applied to 
the numerical results, as indeed is required for presentable results even for 
standard second-order finite difference methods. 

A third application of the SFD method has been to a set of inviscid fluid 
equations with periodic boundary conditions that model some forced gas flows of 
astrophysical interest. The steady-state solutions exhibit strong shocks for 
which the advection steepening is balanced primarily by Coriolis forces rather 
than by pressure forces. The density obtained in one SFD calculation which in- 
cluded a small, explicit viscous term is shown in Figure 3. 
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OUTFLOW BOUNDARY CONDITIONS FOR FLUID DYNAMICS 
Alvin Bayliss* and Eli Turkel 


In many problems in exterior aerodynamics one uses a time-dependent code 
to march the fluid dynamics equations to a steady state. Since these equations 
admit wave-like solutions, a steady state is achieved only because energy radi- 
ates out to infinity. In order to compute the solution in an exterior region 
the domain of integration must be artifically restricted. Mapping the infinite 
region into a finite region is frequently not useful because the waves are not 
resolved near infinity. 

In this work we develop boundary conditions at the artificial surfaces 
which reduce the reflections into the computational domain. In addition, the 
boundary conditions lead to a well posed problem. The conditions are generali- 
zations of those previously der;ived for the wave equation (ref. 1). 

When a subsonic outflow occurs it is necessary to specify one condition 
at the outflow boundary. Frequently the pressure is specified, but this leads 
to large reflections from the boundary. For supersonic flow the top boundary 
is frequently subsonic with respect to the normal velocity. This again requires 
the imposition of one boundary condition. Extrapolating all the variables de- 
grades the accuracy of the solution. 

These new boundary conditions involve time derivatives as well as point 
values of the pressure and velocity fields. Hence, these conditions are easy 
to implement. 

Two test cases were considered. In the first case the Navier-Stokes equa- 
tions were solved with initial conditions that were constant except for a loca- 
lized disturbance. Analytically, the disturbance is convected out of the compu- 
tational domain. This simple problem tests the ability of the boundary condi- 
tion to transmit the disturbance without introducing artificial reflections. 

This problem was run with several different exit boundary conditions. The re- 
sults of Table 1 clearly indicate the improvements obtained. 

In the second case the Navier-Stokes equations with subsonic flow are 
solved over a flat plate. The radiation conditions were imposed at both the 
outflow and top boundaries. The results are presented in Table 2. It is evi- 
dent that the imposition of the correct radiation conditions permits a substan- 
tial reduction in the size of the computational domain as well as in the number 
of iterations required to reach steady state. 

Further studies involving nozzle flow also show substantial improvement 
with the use of the new boundary conditions. 


*ICASE, 505-31-83-01, 804-827-2513 


148 



TABLE 1 


Boundary Condition 


No. Iterations 


Rudy-Strikwerda (ref. 2) 
new 


> 20000 
3383 
1850 


TABLE 2 


b.c. (outflow) 

b.c. (top) 

Position of Top 

No. Iterations 

P = P 

extrapolation 

10 

> 20000 

new 

extrapolation 

10 

13200 

new 

new 

10 

8750 

P = Poo 

extrapolation 

5 

inaccurate 

new 

new 

5 

8800 
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COMPUTATION OF THE FLOW AT A SHUTTLE-TYPE 
WING ELEVON JUNCTURE 

Joanne L. Walsh* 
and 

John C. Strikwerda** 

Some type of sealing arrangement is required to prevent ingress of hot 
boundary layer gas into the unprotected internal wing structure at the 
juncture of the shuttle orbiter wing and eleven. (See left portion of 
figure 1). The cove environment is unknown but potentially hostile, and 
experiments and analyses are underway at Langley to define it. 

In the numerical study a 2D mathematical model using the full 
two-dimensional compressible Navier-Stokes equations for laminar supersonic 
flow at the wing-elevon cove juncture has been developed. An existing 
computer program (Ref. 1) which runs on the Cyber 203 computer was modified to 
do the study. An explicit finite difference technique is used to solve the 
governing equations. The computational domain is discretized into over 6500 
grid points and at each point two velocity vectors, density, and temperature 
are determined. Numerical size of the problem and a small computational time 
step restricted the model to the area shown in the lower right portion of 
figure 1. To accelerate the convergence to steady state, the stability 
condition at each point is used to determine a local time step which is used 
to advance the solution. 

The analysis has defined the complex interaction between the external 
boundary layer flow and the flow in the cove. The basic flow phenomena for a 
nominal leak rate of 3 kg/m^-sec (0.6 Ib/ft^-sec) is illustrated in the lower 
right of the figure. As the flow expands off the wing surface, a 
recirculating eddy develops at the entrance to the cove. A small portion of 
the external boundary layer flow between the body and the dividing streamline 
is ingested through the cove. The remainder of the boundary layer is 
compressed by the deflected eleven surface and results in a shock. The leak 
in the cove is simulated by specifying the pressure at the cove exit to be 
less than the free-stream pressure. The flow description obtained from the 
detailed analysis provides the required insight into the overall phenomena to 
permit simplifying assumptions that make an aerothermostructural analysis of 
the wing-el evon-cove structure possible. The flow analysis is being continued 
to study the effects of different wing surface temperatures and the effects of 
a range of leak sizes, including the no leak condition (i.e. sealed cove), on 
the cove environment. 

References: 

1. Strikwerda, J. C.: A time-split difference scheme for the compressible 

Navier-Stokes equations applied to flows in slotted nozzles. ICASE report 
to appear. 
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SHUTTLE WING-ELEVON JUNCTURE 



Figure 1.- Basic flow description at wing elevon juncture with cove leakage. 



COMPUTATION OF THE FLOW IN SLOTTED NOZZLES 


John C. Strikwerda* 


The steady-state laminar flow in nozzles with suction slots has been com- 
puted using the compressible Navier-Stokes equations. Both two-dimensional and 
axi symmetric flows have been computed. The purpose of these computations was tc 
aid in the design of the quiet wind tunnels being developed at NASA Langley 
Research Center (Ref. 1 and Fig. 1). 

The finite difference scheme is an explicit time-split scheme with three 
splittings. The equations are split into two one-dimensional hyperbolic schemes 
for the inviscid terms and a parabolic scheme for the viscous terms. Each of 
these schemes is second-order accurate in space and time and highly vector! za- 
ble. 


The computational grid was generated numerically using a variant of the 
Thompson method. The computational domain was the union of two rectangles 
joined along a portion of one side. 

Flows were computed for Reynolds numbers as high as 20 million with Mach 
numbers of up to 2 at the exit boundary. The initial flow field was based on 
the one-dimensional inviscid theory of Laval nozzles. The steady-state solu- 
tion was computed by solving the time dependent Navier-Stokes equations until a 
steady state was attained. 

The calculations were performed on the CDC STAR-100 vector processor. The 
computer code is written in SL/1 and utilizes 32-bit arithmetic. 

Reference: 

1. Beckwith, I. E. "Development of a High Reynolds Number Quiet Tunnel for 
Transition Research" AIAA Journal, Vol . 13, No. 3, March 1975, pp. 300-306. 
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Figure 1, Illustration of a Half-section of a Nozzle 
with a Suction Slot. 
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A CLASS OF IMPLICIT DIFFERENCE SCHEMES FOR 
COMPRESSIBLE NAVIER-STOKES EQUATIONS 
Milton E. Rose* 


This research investigates weak-element approximations (ref. 1) to mixed 
initial -boundary value problems for systems in conservation form. Employing 
linear approximations to the solution of the equation ut + fx(u)=0 in local 
computational subdomains, an implicit system of equations for the flux f 
results which is solvable by two-point boundary-value methods; the solution 
u is then explicitly updated by a leap-frog scheme. One class of solution 
techniques is furnished by the Keller Box-scheme. For a scalar equation in 
two-dimensions the method has been verified to be second-order accurate using 
an ADI-type solution method for the flux-equations. 

Two methods for extending this approach to the Navier-Stokes equations 
are being considered. One employs quadratic approximation elements; the 
other adapts the preceding discussion to a certain canonical form of the 
Navier-Stokes equations. 

References : 

1. Rose, Milton E. "Weak-Element Approximations to Elliptic Differential 
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COMPLETE VISCOUS FLOWFIELD SOLUTIONS ABOUT 


A BLUNT PARABOLIC BODY IN A SUPERSONIC STREAM 
K, J. Weilmuenster and R. A. Graves, Jr.* 


The calculation of complete viscous compressible flowfields about blunt 
bodies in a supersonic stream has been a desirable, yet unreachable, goal of 
computational fluid dynamics until recent advances in computer hardware 
provided the necessary tools to attack the problem. There have been many 
partial and "patchwork" solutions to the problem which tend to apply to 
specific conditions and which break down when subsonic flow extends from 
the stagnation point to the aft body comer (sonic-comer flow). Also, 
these specialized techniques cannot handle low Reynolds number flows where 
the viscous region is extensive, and where viscous forces predominate. The 
purpose of the work described here is to demonstrate a general technique 
for computing flowfields about axisymmetric blunt bodies in a supersonic 
stream and to use the technique to investigate blunt-body flows over a 
realistic Mach number/ Reynolds number range. 

Solutions are obtained using the physical coordinate system shown in 
figure 1. A complete account of the solution technique can be found in 
reference 1. Briefly, the flow is divided into two parts, an outer inviscid 
region over which the Euler equations are integrated, and a viscous region 
over which the full Navier-Stokes equations are integrated. The two 
solutions are coupled along the line q = in figure 1, and are marched 
together in time to obtain a solution. Solutions are carried out on the 
CDC-CYBER-203 vector processing computer. 

Solutions for the present sonic- comer body have been obtained over the 
range 2 < M^ < 5 and 500 < Rgj^ < 2000 where Rej^ is the free stream 
Reynolds number based on nose radius. The effects of base flow on total drag 
coefficient was foimd to be significant at low Reynolds numbers. Also, the 
growth rate of the base recirculation region has been investigated along with 
wall pressure and shear stress distributions. In figure 2, is shown a Mach 
number contour plot which is representative of the flowfields computed in 
this study. 


Reference : 

1. Weilmuenster, K. J. and Hamilton, H. H., "A Hybridized Method for 
Computing High Reynolds Number Hypersonic Flow About Blunt Bodies," 
NASA TP 1497, October 1979. 
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VISCOUS COMPRESSIBLE FLOW ABOUT BLUNT BODIES USING A 

NUMERICALLY GENERATED ORTHOGONAL COORDINATE SYSTEM 

Randolph A. Graves, Jr.* 

H. Harris Hamilton, II 


Numerical solutions to the Navier-Stokes equations are 
obtained for blunt axisymmetric entry bodies of arbitrary shape 
in supersonic flow. These equations are solved on a finite- 
difference mesh obtained from a simple numerical technique^ 
which generates orthogonal coordinates between arbitrary 
boundaries. The governing equations are solved in time-dependent 
form using Stetter's improved stability three-step predictor- 
corrector method. Since the overall solution is time— dependent, 
the coordinate system is regenerated at each step as dictated by 
the shock wave movement. This concept was successfully 
demonstrated^ for a large number of inviscid flow calculations. 

For the present application, the metric coefficients were 
obtained numerically using fourth order accurate finite-difference 
relations and proved to be totally reliable for the highly 
stretched mesh used to resolve the thin viscous boundary layer. 

Solutions have been obtained for a range of blunt-body nose 
shapes including concavities. A number of solutions have been 
compared with both experimental and computational results with 
very good to excellent agreement obtained. Figure 1 shows a 
series of nose shapes seen using the present technique and 
Figure 2 gives the computed heat transfer for these configura- 
bions. Overall results indicate that the numerically generated 
coordinate system performed exceptionally well and no problems 
were encountered in the coupling of the numerical coordinate 
generator and the fluid dynamic equations. 

References : 

1. Graves, Jr., R. A.: Application of a Numerical Orthogonal 

Coordinate Generator to Axisymmetric Blunt Bodies, 

NASA TM 80131, Oct. 1979. 

2. Hamilton, H. H. and Graves, Jr., R. A.: Application of a 

Numerically Generated Orthogonal Coordinate System to the 
Solution of Inviscid Axisymmetric Supersonic Flow Over 
Blunt Bodies, NASA TP 1619, Jan. 1980. 


*SSD, 506-51-13, 202-755-3277 
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PREDICTION OF 3-D VISCOUS FLOWS IN THE VICINITY OF 
TRANSONIC WIND-TUNNEL WALL SLOTS 
Douglas L. Dwoyer* 


In recent years, a study has been underway to understand the detailed 
flow in the vicinity of transonic wind-tunnel wall slots. The original 
study was conducted by Dr. Charlie Cooke of ODU and consisted of solving 
the 2-D Navier-Stokes equations in a plane normal to the wall slot. Quali- 
tative comparisons between the computed results and experimental measure- 
ment were good enough to encourage development of a full three-dimensional 
model . 

As a result of the two-dimensional study, a joint effort between the 
author and Dr. Cooke was begun to computationally model a slotted transonic 
wind-tunnel test section. The model is based on the compressible split- 
velocity Navier-Stokes technique in which the total velocity is split into 
its rotational and irrotational parts. This method has the advantage of 
allowing the viscous effects to be included in local regions of the flow 
and having the remainder of the flow governed by the compressible potential 
equation. A schematic of the computational model is shown in the figure. 

It is proposed to include the viscous effects only in the immediate vicinity 
of the slots and treat the flow in the test section and plenum as inviscid. 

The computations will be carried out on the CYBER-203 computer and, 
hence, vector processing compatible numerical algorithms are being used. 

The compressible potential equation will be solved using the ZEBRA II algo- 
rithm (ref. 1) while the viscous equations will be solved by the hopscotch 
algorithm. 

References: 

1. South, J. C. , Jr.; Keller, J. D.; and Hafez, M. M. "Vector Processor 
Algorithms for Transonic Flow Calculations" AIAA Paper No. 79-1457. 

July 1979. 


*STAD, 505-31-13, 804-827-2627 
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PREDICTION OF VISCOUS FLOW IN AERODYNAMIC JUNCTURE REGIONS 
WITH THE PARABOLIZED NAVIER-STOKES EQUATIONS 
Douglas L. Dwoyer* 


Understanding of the details of turbulent flows in aerodynamic juncture 
regions such as wing/fuselage junctions is important in the design of proper 
fillets. Due to the complexity of such flows, however, conventional boundary- 
layer methods cannot be used. For this reason, an investigation was begun by 
Dr. A. J. Baker of CMC, Inc., to use the steady, time-averaged, parabolized 
Navier- Stokes (PNS) equations to predict such flows. In this work the PNS 
equations are used to predict the flow in the viscous layer in the corner 
and are interacted with a three-dimensional inviscid calculation in order to 
account for the global effect of the corner flow. 

In this research, Dr. Baker has formulated the appropriate governing 
equations, developed a suitable turbulence model, and devised an interaction 
procedure. The governing equations are the steady time-averaged compressible 
Navier-Stokes equations with streamwise diffusion terms neglected and with 
the pressure gradient vector split into a streamwise component which is im- 
posed from the inviscid solution and a component which corrects the normal 
pressure gradient which is calculated in the PNS solution. The turbulence 
model is based on a parabolized form of the turbulence kinetic energy equation. 

Shown in the figure is a preliminary calculation showing the cross-flow 
velocity vector at the 12-percent chord station in the junction of two para- 
bolic arc airfoils. These data were plotted from the first viscous/inviscid 
interaction. The beginnings of the rollup of the corner vortices can be seen. 


*STAD, 534-02-13, 827-2627 
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Figure 1.- Parabolic arc juncture region flow 


163 


IMPLICIT METHODS FOR THE NAVIER-STOKES EQUATIONS 

•k 

Stephen F. Wornom 


Implicit factored methods, when applied to the Navier-Stokes equations, 
converge very rapidly as long as the flow is predominantly viscous. The case 
of a shock-wave interacting with a boundary layer on a plate is a good example. 
Results obtained via the method reported in reference 1 show acceptable con- 
vergence in 100 time steps with u^At/Ax = 100 during which the residual error 
is reduced 4 orders of magnitude. On the other hand, for flows that are not 
viscous dominated, the convergence rate is much slower. Consider flow over a 
2-D airfoil using the Euler equations. Results obtained using the method of 
reference 2 are shown in figure 1 where it can be seen that acceptable con- 
vergence has not been reached after 300 time steps with u^At/Ax = 5 in which 
the residual error was reduced one order of magnitude. °° 

The purpose of this research is to develop an implicit two-point method 
for the Euler equations which avoids two weak points of the three-point method 
used in reference 2. The weak points of reference 2 are: (1) The difference 

scheme, when applied to the Euler equations, weakly couples the odd and even 
points and as such tends to produce oscillatory results and (2) the resulting 
difference equations require more boundary conditions than the P.D.E. 's being 
solved. These improvements are part of the overall research goal to develop 
an efficient method for solving the Euler equations and, likewise, the Navier- 
Stokes equations. The related research areas are the multi-grid and flux- 
splitting methods. 

References : 

1. Beam, Richard M. ; and Warming, R. F. "An Implicit Factored Scheme for the 
Compressible Navier-Stokes Equations" AIAA J., vol . 16, no. 4, April 1978, 
pp. 393-402. 

2. Beam, Richard M.; and Warming, R. F. "An Implicit Finite-Difference 
Algorithm for Hyperbolic Systems in Conservation-Law Form" 0. Comp. 

Physics, vol. 22, 2976, pp. 87-110. 


*STAD, 505-31-13, 
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VECTOR COMPUTER APPLICATIONS SOFTWARE 


FOR THE SOLUTION OF THE COMPRESSIBLE 
NAVIER-STOKES EQUATIONS 
Robert E. Smith and Joan I, Pitts 


A three-dimensional Navier-Stokes solver for compressible 
laminar flow based on a proven nimerical technique (MacCormack 
time-split algorithm) has been developed for the CYBER 203 
computer. This explicit technique is used because of its 
robust characteristics and because it conforms to the CYBER 203 
architecture. There is a large data base consisting of five 
state variables and nine variables for geometric description 
at each grid point. The virtual memory of the CYBER 203 requires 
that considerable attention be given to program organization. 

The data base has been interleaved using a four -dimensional 
array and intermediate storage is used to apply the finite 
difference algorithm (fig. 1). 

In developing the Navier-Stokes solver, primary consideration 
has been given to the capability to solve a wide class of 
problems with a minimum of programing changes. This has been 
accomplished by programing the complete transformed equations 
of motion including all nine elements of the Jacobian matrix 
(metric data) . Supplying the metric data from a grid generation 
technique and programing the boundary conditions for a given 
problem (separate subroutine, fig. 2) allows virtually any 
laminar fluid flow problem to be attacked. Experience, however, 
indicates that supersonic flow is most feasibly simulated with 
the MacCormack technique. 

The Navier-Stokes solver is written in the SL/1 programing 
language and uses 32 bit arithmetic. There are approximately 
two million words of incore memory and a computational rate of 

5.5 X 10"^ seconds per grid point per time step is observed. 

This is 26 times faster than the same technique written for 
the CDC 7600. References 2 and 3 present solutions obtained 
with the Navier-Stokes solver. No significant accuracy loss 
has been observed using the 32 bit word size. The program 
is independent of the grid generation technique and the metric 
data from any technique can be used by the program. 
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FINITE ELEMENT TECHNIQUES FOR HIGH 
REYNOLDS NUMBER FLOW PREDICTION 
A. J. Baker* 


A particularly difficult problem class in computational fluid dynamics is 
flowfield prediction at large Reynolds number. The governing Navier-Stokes 
equations are strongly non-linear elliptic boundary value description, but the 
importance of the viscosity term can be completely dwarfed by the non-linear 
convection terms everywhere away from a wall. Furthermore, the continuity 
equation possesses no viscosity-like term, hence uniformly exhibits the hyper- 
bolic description pervading the entire flowfield region. 

The technical goal of this project is to extend the developed finite 
element tensor product numerical algorithm to solution of systems of equations, 
in particular the multi-dimensional Navier-Stokes equations for compressible 
flow at large Reynolds number. The primary requirement is to assess formation 
and solution of the tensor matrix product non-linear algebraic equation 
systems produced by this algorithm, with regards to solution stability^ accuracy 
and speed. A second requirement regards evaluation of a non-iterative algor- 
ithm, with its additional core requirements but reduced number of algorithm 
sweeps. Thirdly, is it plausible to utilize multi-dimensional finite element 
matrices for formation of the non-linear algebraic equation system, while 
employing tensor products only for the Jacobian of the Newton iteration algor- 
ithm. 


The basic algorithm statement employs interpolation and a weighted resid- 
uals formulation to approximate the Navier-Stokes system. Letting L(q) and 
£(q) denote respectively the governing equation and the corresponding boundary 
conditions, where q = (p,pu-j) is a generalized dependent variable, the 

M j 

approximate solution q* = Z {N} {Q> is formulated by requiring the error to 

e=l ® 

be orthogonal to the space of functions {N} used for interpolation on the 
finite element domain Rg, i.e. 


,{N}L(q*) + 


„{N}VL(q*) + X {N}S,(q*) = {0} 

Jr" J 8R 


( 1 ) 


Equation (1) is a system of non-linear ordinary differential equations, and 
the term modified by f introduces a dissipation mechanism for error control 
Employing the trapezoidal integration rule produces a non-linear algebraic 
equation system {F>, the solution of which using Newton iteration is of the 
form 


[J]{6Q} = - {F} 


( 2 ) 
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where {6Q} is the iteration vector and [J] the Jacobian of {F}. 

The emphasis in the first phase of this project, initiated 1 March 1980, 
is to construct regular and tensor product forms for {F>, using both the 
conservative and non-conservative equation forms, and to compare algorithm 
performance. The Jacobian is always formed using matrix tensor products for 
solution economy, and its formulation is invariant for the non-iterative 
algorithm as well. The various matrix forms have been derived and the 
corresponding computer code modifications are underway. Preliminary results 
for the rotating cone test case, corresponding to solution of the continuity 
equation alone, indicates the iterative full matrix {F} solution form, and the 
approximate tensor product {F} form, in essential agreement. The former 
required four iterations/step for convergence to 10"“* using the tensor product 
Jacobian. The complete approximate tensor product form is about a factor of 
four faster on the grid sweeps. The next step is to assess formulations of 
the exact {F}, to improve the solution economy for application to the complete 
Navier-Stokes system. This will primarily involve the approximation to the 
initial -value term, which appears the key to solution accuracy. 
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General Aerodynamic Simulation 
Three-Dimensional, Compressible, Navier-Stokes Equations 

* 

Julius E. Harris 

Substantial progress has been made over the past decade in numerically solving 
the Navier-Stokes equations. .Recent developments in coordinate generation 
procedures, for example, boundary- fitted approaches have allowed solutions for 
systems that could not be previously obtained with conventional Cartesian 
grids. The boundary-fitted approach has also been extended to three- 
dimensional flows for cases where two-dimensional boundary- fitted regions are 
joined by piecewise analytic functions to bound the three-dimensional compu- 
tational domain. Most existing programs treat only perfect-gas, laminar 
flows; however, some progress has been made in implementing turbulence closure 
for pressure driven (boundary- layer like) flows. With the exception of a few 
specialized programs (NASA ARC/ILLIAC IV), almost all software to date has 
been developed for scalar-mode processing on systems such as the CDC-7600. 

The introduction of advanced computer architecture (for example; CYBER-203 
at LaRC) presents the challenging opportunity to develop software that could 
take full advantage of the increased processing speed associated with vector- 
pipeling processing as well as virtual memory storage. Accordingly, an in- 
house program was initiated in the following areas: (1) development of an 

algorithm for processing on the CYBER-203; (2) algorithm structure 
to minimize page-fault penalties associated with virtual memory to maximize 
processing speed; (3) development of a data base management procedure to allow 
maximum use of all computer software/hardware procedures at LaRC; (4) develop- 
ment and/or utilization of a general mesh generation program external to the 
algorithm module; (5) implementation of turbulence closure (algebraic to multi- 
equation systems). To date, goals (1) to (3) have been achieved, progress has 
been made on goal (4) and equation systems have been developed for goal (5). 

An explicit, time-split algorithm was developed and programmed as a stand- 
alone modular package. Modifications occur only in the modules that treat 
boundary and initial conditions and grid generation. Variable dimensions are 
used to minimize computer resource requirements. Page-fault penalties are 
minimized through data base interleaving and vector turn-over as the operators 
sweep through the data base. Flow visualization includes line, contour, 
surface, and color plot capabilities. 

As a first series of test cases the program was applied to supersonic flow 
past sharp and filleted corners. The corner fillets were chosen to be circular 
arcs. Numerical results are presented in Figures 1 and 2. Calculations are 
now being made to determine the effect of discontinuities in the second 
derivativeof wall shape on the vorticity field. Current emphasis is focused 
on developing/incorporating a more general procedure for coordinate generation. 


*HSAD, 505-31-13-03, (804) 827-4546 
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UNSTEADY OPEN CAVITY FLOW/TWO-DIMENSIONAL 
INCOMPRESSIBLE NAVIER-STOKES EQUATIONS 

Joshua C. Anyiwo 

A computer code is being developed to numerically study unsteady, incompressi- 
ble shear flow over a series of open cavities of arbitrary aspect ratios and 
spacing (Fig. 1). 

The computer code processes the vorticity equation (using a simple ADI method), 
and the stream function and pressure Poisson equations (using an SOR method). 
Reference 1. For generality all velocities are non-dimensionalized with u 
and all lengths with A where: ” 

u = (wall shear stress/fluid density)^*^^, at some convenient 
^ reference position; here, the inflow plane; 

A = k {Z A^} {Z (£^ + w^.)}; i = 1, 2, 3 ... 

A^. = aspect ratio of i-th cavity 

= width of i-th cavity 

w^. = spacing between i-th and (I + l)-th cavities 
k = arbitrary numerical constant 

Particular attention will be directed towards the effects of vortex shedding 
and inflow velocity profile on the general evolution of vorticity in the flow 
system. The possibility of applying these results to studies in potential 
drag reduction experiments will be explored. 


Reference: 

1. Numerical Studies of Incompressible Viscous Flow in a Driven Cavity. 
NASA SP-378, 1975. 


*HSAD 505-31-13-03, (804) 827-2806 
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A TWO-DIMENSIONAL NAVIER-STOKES SOLVER USING SPECTRAL METHODS 
FOR FLOWS OVER MOVING WALL GEOMETRIES 

R. Balasubramanian and Steven A. Orszag* 

Thestudy of mechanisms of possible drag reduction effects in flows 
over moving wavy walls as well as the study of mechanisms of generation of 
surface waves by wind requires understanding of the detailed flows over wavy 
surfaces. We have developed a computer code based on spectral methods for 
the study of two-dimensional incompressible flows in wavy geometries. The 
full time-dependent two-dimensional Navier-Stokes equations are solved using 
spectral methods to achieve high spatial accuracy and high-order time-split- 
ting methods together with conformal mapping methods to allow simulation of 
flow over steep waves. 

Applications have been made to laminar flows over wavy surfaces. Some 
experimental results have been given by Kachanov Kozlov, Kotjolkin, Levchenko 
and Rudnitsky at a Reynolds number of 1.1 x 10^. The results obtained using 
our code are in excellent agreement with experiment. 

In addition comparison of the results of the computer codes with the 
theoretical predictions of Benjamin show that the theoretical analysis is 
valid for small values of ka and 6/a, where a is the amplitude of the 
wall displacement, k is its wavenumber, and 6 is the boundary layer 
thickness. 

Our main motivation for the present calculations is to study the pos- 
sibility for drag reduction in flow over non-sinusoidal periodic surfaces. 

The total drag for the surfaces are compared to the total drag over a sinu- 
soidal surface of the same wave length and amplitude. 

We have also studied the laminar flow over a localized bump. The re- 
sults are compared with the results obtained by Taylor and Gent. This study 
should be useful in understanding the effect of local roughness on flows. 

In the area of turbulent flows, the computer codes have been used to 
study (1) channel flows with one wall of the channel replaced by a v/avy 
surface and (2) boundary layer flows over a wavy wall. One of the key 
problems addressed in this study is the kind of turbulence modeling necessary 
for proper simulation of turbulent flow over wavy walls. 

For the boundary layer flow problem, both curvature and pressure gradi- 
ent effects are of importance. Results are obtained using the Beckwith- 
Bushnell model. Current work is being done to improve the turbulence model- 
ing in the boundary layer case. 


*Camb ridge Hydrodynamics, Inc. 
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TIME-DEPENDENT COORDINATE SYSTEMS 
FOR THE NUMERICAL SOLUTION OF FLUID FLOW 

* 

Liviu Lustman 

Powerful algorithms have been used recently to map complicated geometri- 
cal regions over standard squares or cubes. In the corresponding new 
coordinates, uniform grids may be used efficiently, as well as techniques 
necessitating simple geometry, such as spectral methods. However, the grid 
generation is usually effected only once, before starting the solution of the 
differential system. The coordinates thus obtained are appropriate for 
boundary layer resolution or sharp corners, but are essentially static and 
probably not applicable to time-dependent problems. It is a serious challenge 
to design a time-dependent grid, which will map a uniform computational mesh, 
at various times, onto surfaces in the physical domain fitted to wave phe- 
nomena - sharp gradients or shock discontinuities propagating in space. 

Proper bunching of these curvilinear coordinates may solve the shock resolu- 
tion problem, since - as numerical experience shows - Gibbs phenomena and 
shock smearing usually occurs not on a fixed length scale, but on a fixed 
number of mesh points. Clearly, such an ambitious project must overcome 
several difficulties: 

(1) Grid generation is time consuming, if done by the elliptic system method 
of Thompson et al . Simpler geometrical methods are hard to match to 
flow properties in the physical domain. 

(2) After being expressed in the new coordinates, the transformed equations 
will contain new terms due to coordinate evolution, posing stability and 
accuracy problems. 

However, recent results of K. Miller and R. N. Miller show that a very simple 
node algorithm - which is moreover independent of the physical system to be 
solved - is quite efficient in one-dimensional problems. 

I propose to study similarly conceived algorithms for sample two- 
dimensional computations, in order to derive a gradient sensitive, time- 
dependent coordinate generator, and to test its usefulness for various 
numerical procedures. 


^Miller, Keith, Miller, Robert N., "Moving Finite Elements," to appear. 

^Eisemam, Peter,; J.C.P. 26(1978)307, J.C.P. 33(1979)118 

^Thompson, J., Thames, F., Mastin, C.; J.C.P. 15(1974)299, J.C.P. 24 (1977) 
274, NASA-CR 2729. 

Old Dominion University 
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PREDICTION OF THE TRANSONIC FLOW OVER AXISYMMETRIC 
BOATTAIL NOZZLES USING A NAVIER-STOKES CODE DEVELOPED 
FOR THE STAR COMPUTER 
R. C. Swanson* 


The time-dependent Navier-Stokes equations in mass-averaged variables 
are solved for transonic flow over axisymmetric boattail -plume simulator 
configurations. Numerical solution of these equations is accomplished with 
the explicit finite difference algorithm of MacCormack. A grid subcycling 
procedure and computer code vector! zati on are used to improve computational 
efficiency. The two-layer algebraic turbulence models of Cebeci-Smith (C-S) 
and Baldwin-Lomax (B-L) are employed for investigating turbulence closure. 

Two relaxation models based on these baseline models are also considered. 

In this study solutions for the Reubush (Ref. 1) configuration 1 after- 
body, which has a boattail terminal angle of 34 degrees, have been computed. 

For all cases the flow on the boattail is highly separated. In figure (1) 
the variation of the surface pressure coefficient Cp with the nondimensional ized 
axial distance z/D^ at the Mach numbers of 1.3 and ^0.8 is shown. Results 
obtained with each of the four eddy viscosity models are given. The predic- 
tions with the base-line turbulence models show poor agreement with experiment. 
This is not surprising since these models do not account for upstream history 
effects. That is, if the flow is disturbed in some way (i.e., severe adverse 
pressure gradient) the turbulence field does not respond immediately. However, 
it does retain memory of the event. The relaxation models attempt to account 
for this delayed response through a relaxation length scale (X). The Relax 
(C-S) model solutions compare well with the data except for the separated flow 
region. In both the supersonic and subsonic cases the results determined with 
the Relax (B-L) model capture the pressure plateau region. These pressure 
distributions suggests that the relaxation turbulence models have potential 
to accuratley predict such boattail flows. 

Reference: 

Reubush, D. E., "Experimental Study of the Effectiveness of Cylindrical 
Plume Simulators for Predicting Jet-On Boattail Drag at Mach Numbers Up 
to 1.30," NASA TN D-7795, November 1974. 


*HSAD,. 505-32-13, 804-827-2673 
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configuration 1. 



SPLINE COLLOCATION FOR VISCOUS FLOW OVER 
AFTERBODY CONFIGURATIONS 
S. J. Rubin* and R. C. Swanson** 


For a number of years solutions for inviscid transonic flow have been 
obtained efficiently by solving a potential flow equation with some type of 
relaxation procedure. Also, the viscous regions of such flows have been 
determined with efficient solvers of the boundary-layer equations. In this 
work some of the principles of these techniques are combined to solve the 
Navier-Stokes equations. 

The subsonic compressible Navier-Stokes equations are solved using a 
velocity field decomposition procedure. That is, the streamline velocity 
component is represented by a sum of potential and viscous flow contributions. 

In this way the boundary-layer and potential flow behaviors are recovered in 
their respective regions of validity. The fully coupled flow equations, which 
are written in a body-fitted coordinate system, are solved with Stone's 
strongly implicit method. The time derivative merely provides an artiface for 
obtaining the steady state solution. In incompressible computations, both 
unseparated and separated, with this method (Ref. 1) sufficiently large time 
steps have been achieved to collapse this procedure to a relaxation process. 

This characteristic is an objective in the compressible calculations of this 
effort. To further enhance computational efficiency or to increase numerical 
accuracy a higher-order spline interpolation deferred corrector is incorporated. 

At this point compressible attached flow solutions have been determined 
for a boattail like configuration. Work is in progress to compute separated 
flows over such geometries. 

References : 

Rubin, S. G. and Khosla, P. K. , "Navier-Stokes Calculations With a Coupled 
Strongly Implicit Method - Part 1: Finite-Difference Solutions," AIAA 

Paper 79-0011, January 1979. 


*University of Cincinnati 

**HSAD, 505-32-13. 804-827-2673 
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A TWO-STREAM SOLUTION OF THE NAVIER-STOKES EQUATIONS 
FOR UNINSTALLED NOZZLES 

Michael C. Cline* and Richard G. Wilmoth** 


The need for a fully viscous solution technique for the prediction of 
nozzle afterbody-exhaust jet interactions at subsonic and transonic speeds 
has been amply demonstrated in recent years. The purpose of this work is to 
develop a two-stream solution technique for the two-dimensional, compressible 
Navier-Stokes equations which can be applied to a wide variety of uninstalled 
nozzle flow problems. This technique will provide not only a basic predictive 
tool but will also be used to guide the development of flow-field models for 
the more approximate patched predictions. 

A two-stream, Navier-Stokes computer code (VNAP2) has been developed as 
an extension of a code developed mainly for treating nozzle internal flows 
(ref. 1). The time-dependent flow equations are marched to steady state in 
the computational domain shown in figure 1 using the explicit, two-step 
MacCormack algorithm. Variable grid spacing in both the physical x- and y- 
directions is used to resolve boundary layers and shear layers. The principal 
features of the computer program are listed in figure 2. A special quick- 
solver technique is also provided to reduce the computational time needed in 
regions of small y-grid spacing (such as near walls and shear layers). 

The program is currently being applied to the calculation of subsonic 
external flow-supersonic exhaust jet problems for which an extensive experi- 
mental data base exists. 


References : 

1. Cline, Michael C. : VNAP: A Computer Program for Computation of Two- 

Dimensional, Time-Dependent Compressible, Viscous, Internal Flow. 
Los Alamos Scientific Laboratory Rept. No. LA-7326, 1978. 


* Los Alamos Scientific Laboratory 

** HSAD, 505-32-13, 804-827-2675 
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Figure 1.- Physical and computational flow spaces. 


Figure 2.- Principal features of viscous flow computer 
program. 



A 3-D SOLUTION OF THE NAVIER-STOKES EQUATIONS FOR THE 
FLOW CHARACTERISTICS AND PERFORMANCE OF NONAX I SYMMETRIC 

NOZZLES 

P. D. Thomas* and Lawrence E. Putnam** 

Recent experimental studies have shown that the nonaxi symmetric 
nozzle properly integrated into twin-engine fighter-aircraft configurations 
has the potential for reducing drag, increasing instantaneous maneuverability 
and subsonic agility, providing STOL capability, and reducing fabrication 
and maintenance costs (ref. 1). To date, the design and integration of 
nonaxi symmetric nozzles have been based, primarily, on semi -empirical tech- 
niques and experience. A solution technique for the 3-D Navier-Stokes 
equations has, therefore, been developed (ref. 2) to predict the performance 
and characteristics of the viscous flow field for uninstalled 3-D nozzle 
configurations at subsonic, transonic, and supersonic speeds. 

The 3-D time-dependent Navier-Stokes equations are solved using the 
Beam and Warming alternating-direction implicit numerical technique. Boundary 
conditions are computed with an implicit technique compatible with the ADI 
technique employed at interior points of the flow region. The equations 
are formulated and solved in a boundary-conforming curvilinear coordinate 
system. The curvilinear coordinate system and computational grid is generated 
numerically as the solution to an elliptic boundary valve problem. A new 
method has been developed that automatically adjusts the elliptic system 
so that the interior grid spacing is controlled directly by the user defined 
grid spacing on the boundary of the flow region. The viscous terms are 
modeled using a modification of the Baldwin and Lomax method. 

A comparison of predicted and experimental surface pressures on the 
center line of a flap of a convergent-divergent rectangular nozzle operating 
at a nozzle pressure ratio of about 3 is shown in figure 1. For these cal- 
culations, the flow was assumed to be two-dimensional. The computational 
grid used is shown in figure 2 and a density contour plot is presented as 
figure 3. The predicted surface pressure are in very good agreement with the 
experimental data. The theory smears the shock located near the throat 
of the nozzle with a resultant decrease in prediction accuracy. Further 
studies are underway to determine the capabilities of the solution technique 
for predicting the flow characteristics of 3-D nozzle configurations. 

References: 

1. Capone, Francis J.: The Nonaxi symmetric Nozzle-It is For Real. AIAA 

Paper No. 79-1810, Aug. 1979. 

2. Thomas, P. D.: Numerical Method For Predicting Flow Characteristics and 

Performance of Nonaxi symmetric Nozzles-Theory. NASA CR-3147, 1979. 
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BOUNDARY CONDITIONS FOR SUBSONIC COMPRESSIBLE 
NAVIER-STOKES CALCULATIONS 
David H. Rudy* and John C. Strikwerda** 


A systematic study has been undertaken of inflow and outflow boundary 
conditions for the numerical solution of the compressible Navier-Stokes 
equations. Combinations of representative inflow and outflow boundary 
conditions were applied in the solution of the subsonic flow over a flat 
plate in a finite computational domain. These boundary conditions were 
evaluated in terms of their effect on the accuracy of the solution and 
the rate of convergence to a steady state. 

An outgrowth of this study was the development of a non-reflecting 
outflow boundary condition which has been shown to be effective in reducing 
reflections from subsonic outflow boundaries. These reflections can have a 
large effect on the rate of convergence of the solution to a steady state, 
and substantial improvements in the rate of convergence to steady state were 
obtained using this non-reflecting boundary condition. In computations of 
transient solutions this new boundary condition is useful to reduce the 
magnitude of the non-physical reflections from the outflow boundary. 

The non-reflecting boundary condition involves a parameter whose 
optimal value is estimated using an analysis of a simplified set of equations 

References : 

1. Rudy, D. and Strikwerda, J. "A Non-Reflecting Outflow Boundary Condition 
for Subsonic Navier-Stokes Calculations" J. Comp. Phys., to appear. 

2. Rudy, D. and Strikwerda, J. "A Non-Reflecting Outflow Boundary Condition 
for Subsonic Navier-Stokes Calculations" ICASE Report 79-2, January, 1979. 

3. Rudy, D. and Strikwerda, J. "Boundary Conditions for Subsonic Compressibl 
Navier-Stokes Calculations" ICASE Report 79-18, August, 1979. 
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NUMERICAL SIMULATION OF HIGH REYNOLDS NUMBER. 
THREE DIMENSIONAL, COMPRESSIBLE FLOW 
M. Y. Hussaini* and Steven A. Orszag 


We have developed a high resolution mixed finite-difference-spectral code 
for the direct numerical simulation of high Reynolds number three-dimensional 
compressible flow. The code runs with up to 64 x 32 x 32 degrees of freedom on 
the STAR-100 computer in about 5 sec. per time step. Semi -implicit time split- 
ting methods are used to achieve second-order accuracy in time while maintaining 
stability on convective time scales for low Mach number runs. 

The code is being applied to the study of the inviscid stability of boun- 
dary and free-shear layers. The code is validated by comparison with known 
linear results of Mack and others. Finite amplitude nonlinear effects in plane- 
parallel compressible shear flows are also studied. The role of three-dimen- 
sionality is examined as a function of Mach number in order to compare mecha- 
nisms of transition to turbulence in compressible and incompressible flows. 
Detailed flow properties are being studied. Finally, the decay of three- 
dimensional homogeneous, compressible turbulence is being studied. Effects of 
compressibility on turbulence spectra are being investigated. 


*ICASE, 505-31-83-01, 804-827-2513 
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Parabolic Approximation 


Methods in which the steady N-S equations are approximated in such a manner 
that a predominant spatial "marching" direction is used. Not including 
Boundary-Layer methods. 
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APPLICATION OF PARABOLIC ANALYSIS TO 
PRACTICAL SCRAMJET FLOW FIELDS 
John S. Evans* 


A vital part of the effort to develop a supersonic combustion ramjet 
engine (scramjet) is the capability for theoretical calculations of flow 
properties in the combustor, where hydrogen is injected into and burned in 
a supersonic airstream. In various parts of the engine the flow includes 
shock waves, expansions, turbulent mixing of gas streams, complex boundaries, 
recirculating flows, and chemical reactions. No computer code exists which 
can handle such an array of difficult problems simultaneously, and, even if 
there were such a program available, it would be too expensive and cumbersome 
for general use. 

There are many situations in supersonic combustion ramjets in which 
simple analyses will suffice. If the flow has a predominant flow direction, 
and if it is supersonic, the parabolized Navier-Stokes equations adequately 
characterize the flow field. These conditions are often satisfied in the 
downstream portion of a scramjet combustor, where the near-field effects of 
fuel injection have disappeared. 

Calculated profiles (ref. 1) typical of such a downstream combustor flow 
field are compared in figure 1 with experimental results. Helium was injected 
from a fuel injector strut parallel to the supersonic airstream, and profiles 
were taken at the center exit plane of helium mass fraction, temperature, and 
flow velocity. The agreement between the data and the calculated profiles is 
fairly good and is representative of the kind of predictive capability needed 
for scramjet combustor development. 

References : 
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Modeling for Scramjet Engines." Paper presented at the AGARD 54th (B) 
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COMBUSTOR ANALYSIS CAPABILITY 



•3-D TECHNIQUE FOR DOWNSTREAM COMBUSTOR FLOW HAS BEEN DEVELOPED 
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Figure 1. Simulated fuel injection profiles (helium); theoretical and experimental 


TURBULENT MIXING AND REACTION IN THREE-DIMENSIONAL PARTIALLY 
ELLIPTIC FLOWS OF SUPERSONIC COMBUSTORS 
R. Clayton Rogers* 


The ability to computationally predict details of the turbulent, chemi- 
cally reactive flows in supersonic combustor configurations has been steadily 
improving over the past ten years. Until recently, these flows have been 
modeled using parabolic assumptions so that the computer code could employ a 
spatially marching solution algorithm. As long as the combustor flow is super 
sonic with a predominant flow direction, at conditions such that chemical 
reaction rates are fast, and the turbulence is properly initialized, these 
codes can be used successfully. However, in certain regions of the combustor, 
such as near fuel injectors where flow separation and recirculation may occur 
and in zones where chemical reactions may produce local regions of embedded 
subsonic flow, the flow may be elliptic. Even with the more powerful computer 
hardware now available, and improved numerical algorithms, the solution of the 
fully elliptic form of the governing equations is intractable at the present 
time. Therefore, a partially elliptic scheme has been developed to handle 
flows which are elliptic due to local subsonic flow but do not have separation 
and recirculation. 

The present computer code was adapted from the SHIP computer program 
(ref. 1). SHIP is capable of computing the three-dimensional parabolic flow 
of the turbulent mixing and combustion of hydrogen in a supersonic airstream. 
In the partially elliptic version, all variables except static-pressure are 
computed as if the flow was purely parabolic. The static pressure, however, 
is stored in a three-dimensional array and is used in subsequent iterations to 
determine the approximate pressure gradient source term in the momentum equa- 
tions. The numerical procedure combines the marching integration scheme of 
SHIP with an iterative scheme to solve the elliptic pressure field. Thus, the 
procedure is termed "partially" elliptic. Additional details of this code 
may be found in reference 2. 

Current plans are to evaluate this partially elliptic code with super- 
sonic combustion data in which local subsonic flow regions are expected. Such 
evaluations will lead to further development of the iterative procedure and a 
possible criterion for convergence of the iterative solution. 

References : 
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BLUNT-BODY VISCOUS-SHOCK- LAYER ANALYSIS THAT INCLUDES 


TURBULENCE, MASS INJECTION, AND RADIATION 
James N. Moss* 


A viscous-shock- layer analysis has been developed in-house at the 
Langley Research Center to specifically study the flow field and heating 
environment experienced by planetary probes. The computer program for this 
analysis is called HYVIS. The motivation for developing the program was the 
fact that the aerotheimal environment encountered by entry probes cannot be 
duplicated in existing ground experimental facilities. Consequently, probe 
heat shield designs must rely extensively on numerical predictions. 

Because of the complexity of the flow about entry probes, the numerical 
solutions must account for physical phenomena such as ablation injection, 
radiation transport, reacting chemistry, and turbulence. 

The HYVIS program has been written in CDC FTN Fortran Extended and has 
been designed for optimum execution in the CDC CYBER environment. The flow 
equations are obtained from the steady-state Navier-Stokes equations by 
retaining terms up to second-order in the inverse square root of a 
Reynolds number (Ref. 1) . Consequently, one set of equations are solved 
that are uniformly valid throughout the shock layer. Real-gas effects 
are accounted for by treating the gas as an equilibrium gas mixture where 
the equilibrium composition is determined by a free energy minimization 
procedure. Radiation transport is calculated with a nongray radiation model 
that accounts for molecular band, atomic line, and continuum transitions. 

Also, the flow may be either laminar or turbulent. Turbulence is calculated 
by using a two- layer eddy-viscosity model. For solutions with mass 
injection, the injection rate may be either specified or calculated for 
quasi-steady ablation. 

Figure 1 presents results of calculations for a probe (35" hyperboloid) 
entering the atmosphere of Jupiter. The effects of both turbulence and 
coupled (the mass injection rate m is governed by the net surface heating) 
ablation injection on radiative [Fig. 1(a)] and convective [Fig. 2(b)] heating 
are demonstrated. Results such as these are used to establish correlations 
that can be used in approximate engineering calculations. 

Reference: 

1. Moss, J. J., "A Study of the Aerothermal Entry Environment for the 
Galileo Probe," AIAA Paper 79-1081, June 1979. 
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(a) Radiative heating. 

Fig, 1 Effect of coupled carbon-phenolic injection and 
turbulence on surface heating for Jupiter entry 
(freestream conditions: velocity = 39,29 km/s, 

density = 4,36 ^ 10“^ kg/m^, and Mach number 
= 47,6, The nose radius Rjq = 0,3112 m). 



(b) Convective heating 
Fig, 1 Concluded, 


197 


PRELIMINARY THERMAL ANALYSIS FOR SATURN ENTRY 


E. V. Zoby and J. N. Moss* 


The planetary exploration program of NASA includes possible missions 
to the planet Saturn. Since experimental facilities cannot adequately 
duplicate the aerodynamic and thermal conditions expected during Saturn 
entry, much of the required information must be obtained from computational 
analyses. The purposes of this investigation are to scope the thermal 
problem for Saturn entry as well as assess the adequacy of current 
numerical techniques. For a realistic evaluation, the analysis is based 
on a recent mission study referred to as the Saturn Orbiter Dual Probe 
(S02P) mission. 

The preliminary thermal study, ^ using viscous-shock- layer (VSL) and 
engineering codes developed for Project Galileo, investigated nonequilibrium 
chemistry effects on the Saturn thermal environment, defined the primary 
heat-transfer mode for heat-shield design, delineated some problem areas 
for future thermal studies, and validated an engineering code^ for parametric 
or design studies. 

The effect of non equilibrium chemistry appears to significantly 
influence only the radiative fluxes with effects localized to the stagnation 
region. This latter result is demonstrated in Figure 1. The freestream 
conditions are based on the nominal Saturn entry (S02P) . Note the radiative 
heating rates based on an equilibrium VSL code^ are approximately 2 MW/m^ 
where the corresponding value for Project Galileo is approximately 500 MW/m^. 
The heat-transfer mode pertinent to the overall heat-shield design is 
convection as demonstrated by the results presented in Table 1. These results 
are given in terms of the integral of the heating rate over the body surface 
area. 
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FIGURE 1. NONEQUILIBRIUM AND EQUILIBRIUM RADIATIVE HEAT-TRANSFER 
DISTRIBUTIONS. 


TABLE 1. EQUILIBRIUM RADIATIVE AND CONVECTIVE POWER^ 


TIME, sec 

/q,dA 


13.00 

0.04 

6.4 

15.00 

0.10 

10.3 

16.75 

0.16 

15.1 

17.75 

0.16 

17.2 

19.00 

0.12 

20.9 


^35° HYPERBOLOID; POWER = MW, m = 0 
^LAMINAR FLOW, WALL TEMPERATURE = 4000 K 
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Turbulence Effects and Turbulence Modeling 


Studies of the effects of using various turbulence models in essentially 
existing Navier-Stokes or Boundary-Layer codes, or methods for modeling 
turbulence, and turbulence simulation. 
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TURBULENCE MODELS FOR SUPERSONIC COMBUSTORS 


P. T. Harsha* and R. B. Edelman* 


Severe and conflicting demands are made upon a computer code for engi- 
neering calculations of fluid-flow properties in supersonic combustors. In 
each step solutions must be found for several partial differential equations, 
grid-spacing must be fine enough for the resolution of gradients, and the step 
size must be small enough to avoid inaccurate or diverging answers. There is 
a motivation to make all parts of the code as simple and efficient as possible 
so as to minimize computer time and cost. Typical of the kinds of compromises 
which must be made is the selection of a turbulence model. 

The transport of energy, momentum, and mass in the turbulent flow is far 
too complex to calculate in detail and must be modeled. Mixing length models 
(ref. 1) of varying degrees of complexity serve well for boundary layers near 
walls. For free shear layers with a dominant flow direction, the 2-equation 
k-e model (ref. 2) is popular because the dissipation length scale is automati- 
cally calculated and because it models turbulence properties fairly accurately 
for a wide range of flows without changing the empirical constants. Beyond 
these cases, there is no consensus as to the best model, but there is a need 
to determine suitable models to use in complex geometry, in recirculating flow, 
and in reacting flow. 

An investigation is being carried out under Contract NASl-15988 to de- 
termine the best turbulence model (or models) for use in scramjet combustor 
calculations. Present plans call for tests of multiple dissipation length 
scale k-e models (ref. 3), use of an algebraic Reynolds stress model (ref. 4), 
combinations of both of these, and, finally, use of a full Reynolds stress 
model. The models will be tested with 2-D and 3-D parabolic and elliptic 
programs. To the extent that data are available, the models will be adjusted 
so as to maximize agreement with data. To be useful for engineering calcula- 
tions the models chosen cannot be so complex as to increase computation costs 
by an unreasonable amount. 
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MERGING OF TOROIDAL VORTICES 


C. H. Liu* 


Vortex rings have been of interest for many years due to their widespread 
appearance in axi symmetric flows. Many fluid dynamic phenomena can be 
explained using models consisting of vortex rings moving in an inviscid stream. 
Asymptotic solutions for the interaction of free vortices with finite viscous 
vortical cores were based on the assumption that the ratio of the effective 
viscous core size over a certain length scale is small. It is well known 
that the asymptotic solutions may remain as good approximations even when 
this ratio is of the order of unity. In order to establish the accuracy of 
the asymptotic solutions, finite difference approximations to the axisymetric, 
unsteady incompressible Navier-Stokes equations are constructed to study 
the merging of vortex rings. 

An initial distribution of vorticity, representing one or more vortex 
rings, is specified and then allowed to evolve by numerical solution of the 
governing equations. Since the vorticity decays exponentially as the square 
of the distance from a vortex, the vorticity can be assigned zero on the 
boundary which is at sufficient distance from the vortex. The coordinates 
are moving with instantaneous mean velocity of the vortex rings relative to 
the inertial coordinates so that the free vortices will remain inside an 
assigned boundary. 

Numerical results are obtained for the self-merging of a single vortex 
ring which happens as the inner' radius of the ring shrinks to zero. Figure 1 
compares the asymptotic solutions for the decay of the maximum vorticity 
with the numerical results for Re = 12.56. It is clear that the asymptotic 
solution remains quite accurate even at t - 16. The error is caused by 
the nonlinear terms, therefore, the disagreement with the numerical results 
will increase as Re increases. 
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FINITE-ELEMENT SOLUTION ALGORITHM FOR AEROACOUSTIC 
TRAILING EDGE FLOW PREDICTION 
A. 0. Baker 


Turbulent boundary layer flows departing the trailing edge of an aero- 
dynamic surface are experimentally verified to be strong sources of noise. 
Alteration of the aerodynamic surface can alter far-field intensity by local 
absorption as well as modification to the turbulent flow structure prior to 
departing the trailing edge. A considerable volume of experimental data has 
been obtained for the case of a jet flow directed over a surface. Attention is 
now turning to aeroacoustic analysis of the near-field wake associated with 
airfoil flowfields induced solely by forward flight. 

A theoretical analysis is presented yielding sets of partial differential 
equations for determination of turbulent aerodynamic flowfields in the vicinity 
of an airfoil trailing edge. A four-phase interaction algorithm is derived to 
facilitate the analysis. Following input, the first computational phase is 
an elementary viscous-corrected two-dimensional potential flow solution 
yielding an estimate of the inviscid-flow induced pressure distribution. The 
next phase involves solution of the turbulent two-dimensional boundary layer 
equations up to the trailing edge, with transition to a two-dimensional 
parabolic Navier-Stokes equation system describing the near-wake merging of the 
upper and lower surface boundary layers. An iteration provides refinement of 
the potential flow induced pressure coupling to the viscous flow solutions. 

The final phase, if desired, is a complete two-dimensional Navier-Stokes 
analysis of the wake flow in the immediate vicinity of the trailing edge. 

A comparison of computed mean velocity, U^, with experimental data from 
reference 1 downstream of the trailing edge for the airfoil NACA 63-012 at 
zero angle of attack is shown in figure 1. The computed results were obtained 
using the two-dimensional parabolic Navier-Stokes formulation and compare quite 
favorably with data except in the region where the upper and lower streams 
merge. 

This work has been conducted by Computational Mechanics Consultants under 
NASA contract NAS1-14855/M0D. 1. 

Reference: 
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ANRD, 505-10-43-01, 804-827-2617 


206 







5 T M 

%c 

Oj 

i70057 


10ZZ2 

■f" 

1.0'130 

X 

1. 0650 

ilJ 

1.0910 

0 

1. 1000 

ft « • • • 

DATA 


Transverse Coordinate 


ig. 1 


U-, Mean Velocity Profile 



TIME DEPENDENT NUMERICAL SOLUTION OF TURBULENT JET FLOWS 


USING REYNOLDS STRESS CLOSURE MODELS 
T. B. Gatski 


In the study of turbulent flows, it has become common practice to numeri- 
cally solve a set of conservation equations for stationary flow variables to 
determine the mean characteristics of each flow. The level of sophistication 
utilized in the closure of the governing differential equations is determined 
by which flow variables needed to be known and the flow geometry. One can 
classify these closure approximations according to the number of differential 
equations used to describe the required Reynolds stresses. 

In the study of sound generated by free turbulent shear flows, it is 
becoming apparent that the large coherent scales of motion in the flow play a 
significant role in the sound generation process, either directly or through 
an interaction with the smaller component scales of motion. This fact has led 
to experimental investigations which have attempted to enhance these coherent 
scales by acoustic excitation in order to better understand their specific 
role. If the coherent motion mentioned earlier can be modeled in a determin- 
istic manner then a time history of certain scales of motion can be computed. 

In this vein, consider a triple decomposition of the flow field into a 
time mean component, a phase-averaged component, and a small scale random 
component. Phase-averaged equations for the vorticity, stream function and 
(phase-averaged) turbulent stresses are formulated and the entire system 
solved numerically in a radially stretched cylindrical coordinate system from 
the jet exit to eight diameters downstream. The transport equations for the 
phase-averaged stresses are those obtained by Launder, et al . (ref. 1) for the 
time-averaged Reynolds stresses with the constants adjusted for the forced 
flow under consideration. 

The numerical experiment under consideration is the modeling of a forced 
axisymmetric turbulent jet. The initial solution stage consists of the 
solution of the unforced jet problem to provide a consistent set of initial 
conditions for the forced case. Since the governing differential equations 
are time dependent, the initial solution phase requires that the equations be 
computed through to a steady-state solution. This allows for validation of 
the differential model and numerical algorithm as well as insight into the 
transient characteristics of the flow. 

Reference: 
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COMPUTED EVOLUTION OF LARGE SCALE WAVE-LIKE STRUCTURES 
IN TURBULENT JET FLOWS AND THEIR SOUND PRODUCTION 
P. J. Morris 


It is well known in the study of free turbulent jet flows that excitation 
of the jet by a pure. tone causes broadband jet noise amplification. In order 
to explain this phenomenon, it is necessary to understand the effects of the 
pure tone excitation on the structure of the turbulent flow. Fortunately, 
recent studies of turbulent flows indicate that there is a high degree of 
organization at the larger scales, and, in the case of a free turbulent jet, 
that the organized structure may be adaptable to a wave representation. 

One method of isolating the large scales is to decompose the flow field 
into a time average component, a phase averaged component, and a small scale 
random component. From this triple decomposition, energy equations for the 
respective component scales can tie obtained. The closure problem present in 
such a formulation is handled by way of gradient diffusion hypotheses and 
dimensional analysis considerations in the turbulent energy equation and by 
way of shape assumptions, invariant with respect to downstream distance, in the 
large scale energy equation. 

The time mean continuity and momentum equations along with the mean 
turbulent kinetic energy equation are rewritten in terms of a modified stream 
function coordinate and the system solved using the method of cubic splines. 

The large scale energy equation is solved by taking the large scales as 
instability waves of fixed real frequency whose wavelength and growth rate 
are the eigenvalues of a local viscous stability calculation and integrating 
the energy equation in the radial direction. This reduces the partial dif- 
ferential equation for large scale kinetic energy to an ordinary differential 
equation, in the axial direction, for the magnitude (squared) of the wave 
amplitude. A calculation of the local growth rate, -a-j , as a function of jet 
half width, b, is shown in figure 1 for the potential core region of the jet. 

As is seen, the agreement between the viscous and inviscid solution is very 
good, even in the damped region, aj > 0. The far-field pressure distribution 
produced by the large scale structure can be readily calculated using a 
matched asymptotic expansion method developed by Morris and Tam (ref. 1). 

This work is being conducted at the Pennsylvania State University under 
NASA grant NSG 1580. 

Reference: 

1. Morris, P. J. and Tam, C. K. W. "Near and Far Field Noise from Large-Scale 

Instabilities of Axisymmetric Jets," AIAA 4th Aeroacoustics Conference, 

Atlanta, Georgia, 1977, Paper No. 77-1351. 
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8. Viscous Flow Stability and Transition 


Methods whose primary focus is the stability of viscous flows. 



CALCULATION OF 3-D LAMINAR BOUNDARY 


LAYER TRANSITION IN THE PRESENCE 
OF FREE -STREAM TURBULENCE 
E. Clay Anderson* 


The design suction rate distributions required to maintain laminar flow on 
the supercritical airfoil to be used in the Langley LFC swept— wing experiment 
have been determined using a linearized stability analysis, which does not 
include provisions to account directly for the turbulence level in the appro- 
aching freestream. Recent developments in "higher-order" or two-equation tur- 
bulence models have led to predictive methods that directly account for the 
intensity of turbulence in the freestream. Some of the two— equation turbulence 
models have been applied successfully to predict the influence of freestream 
turbulence levels upon laminar-turbulent transition in two-dimensional flow 
fields. More recently, Dr. Wilcox of DCW Industries has developed a three- 
dimensional computer code incorporating an advanced two-equation turbulence 
model which appears ideally suited for the analysis of the boundary layers on 
the LFC airfoil. Dr. Wilcox has recently demonstrated an ability to perform 
calculations which showed the importance of including the influence of free- 
stream turbulence in laminar flow applications. 

Because of the extreme sensitivity of laminar boundary layers to the free- 
stream turbulence intensity, and the high degree of imprecision involved in all 
current prediction methods, an independent analysis of the LFC-airfoil boundary 
layers in the presence of freestream turbulence levels that are representative 
of the known freestream turbulence levels in the Langley 8'-TPT will be done by 
DCW Industries for NASA Langley. 


*STAD, 534-01-13, 804-827-2627 
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NONLINEAR STABILITY THEORY AND OPTIMUM NUMERICS 
FOR LAMINAR FLOW CONTROL ANALYSIS 


•k 

Steven A. Orszag 


A number of extensions and generalizations of 
the SALLY computer code for stability analysis are 
being studied. First, the extension of SALLY to 
study compressible flows is underway. Second, the 
generalization of SALLY to include nonparallel and 
nonlinear effects is underway. Finally, some new 
efficient techniques for the resolution of critical 
layers are being investigated to further improve 
the efficiency of SALLY. 

SALLY is based on the solution of the Orr- 
Sommerfeld equation using expansions in Chebyshev 
polynomials. For incompressible flow, the critical 
layer lies close to the wing so these expansions 
provide good accuracy with modest resolution. On 
the other hand, in compressible flows, the critical 
layer moves away from the wing so additional resolution 
is required to achieve acceptable accuracy. However, 
we have developed a new method to resolve the critical 
layer which promises significant efficiency improvements 
in both compressible and incompressible environments. 

Nonlinear and nonparallel effects are included 
using a spectral transform method. This technique 
permits efficient implementation of these terms with 
minimal technical complications. 


* 

Cambridge Hydrodynamics, Inc. 





THREE-DIMENSIONAL, COMPRESSIBLE, NONPARALLEL, LINEAR STABILITY THEORY 

FOR LAMINAR FLOW CONTROL ANALYSIS 

Ali H. Nayfeh* 

A theory is developed for the linear stability of three-dimensional 
growing incompressible or compressible boundary layers. The method of 
multiple scales is used to derive partial -differential equations describing 
the temporal and spatial evolucion of the complex amplitudes and wave- 
numbers of the disturbances. In general, these equations are elliptic unless 
certain conditions are satisfied. For a monochromatic disturbance, these 
conditions demand that the ratio of the components of the complex group 
velocity be real, thereby relating the direction of growth of the disturbance 
to the disturbance wave angle. For a non-growing boundary layer, this con- 
dition reduces to da/dB being real, where a and B are the complex wave- 
numbers in the- streamwise and crosswise directions, in agreement with the 
result obtained by using the saddle-point method. For a wavepacket, these 
conditions demand that the components of the complex group velocity be real. 
In all cases, the evolution equations are reduced to inhomogeneous ordinary- 
differential equations along real group velocity directions. 

The theory is being used in conjunction with the e*^ method to develop 
a prediction method for transition. The method will be used to evaluate the 
effects of the proposed suction distributions on laminar flow control over 
swept back supercritical wings. 


Engineering Science and Mechanics Department 
Virginia Polytechnic Institute and State University 
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TWO - AND THREE - DIMENSIONAL BOUNDARY -LAYER 
STABILITY ANALYSIS 


Mujeeb R. Malik and Steven A. Orszag 


* 


We compare several methods of transition 
prediction by linear stability analysis. We use 
the incompressible linear stability code SALLY 
in various ways to study flows over laminar flow 
control wings. In particular, we compare the so- 
called envelope method with wave-packet methods 
to predict transition. We conclude that the 
envelope method is at least as reliable as the 
more complicated wave-packet method. 

The code SALLY is based on numerical 
solution of the Orr-Sommerfeld equation for three- 
dimensional disturbances in boundary layer flows 
using expansions in series of Chebyshev polynomials. 

In its original form, SALLY allowed computation 
of (i) maximum amplification among all wave lengths 
and propagation angles; (ii) amplification at 
fixed wave-length and propagation angles; (iii) 
amplification at fixed frequency and fixed wave 
length; (iv) amplification at fixed wave frequency 
and fixed propagation angle. We have extended the 
SALLY code to calculate the propagation of wave 
packets . 

In the Figure, we compare various transition 
prediction methods applied to data taken at Cranfield 
for a large, untapered, 45° swept half wing mounted 
as a dorsal fin upon the mid-upper fuselage of an 
Avro Lancaster airplane. In this case, the angle of 
attack is 0° and experimental transition occurred at 
x/c = 7%. The envelope method gives an N factor 
of 10.8 at a frequency of 1000 Hz. The wave packet 
method gives a maximum N factor of 10.5 at a 
frequency of 1200 Hz, which is close to the prediction 
of the envelope method. A conservative wave packet 
approximation gives an N factor of 8.6 rather than 
10 . 5 . 


Systems and Applied Sceinces , Inc. and Cambridge Hydrodynamics, 
Inc. 
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A plot of N versus precent of chord x/c for various methods applied to a 45 
swept wing mounted on an Avro Lancaster airplane at 0° angle of attack. 
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ON THE STABILITY OF THE FREE SHEAR LAYER FOR AN 
AXIAL SYMMETRIC JET AND PLUG FLOW JETS 

Steven A. Orszag*, Lucio Maestrello**, and 
Stan Lamkin*** 


The stability characteristics of inviscid axial symmetric jets and plug 
jets are being investigated to explain recent experimental observations by 
Maestrello and collaborators on the effects of plugs in the reduction of 
turbulent jet noise production. Both compressible and incompressible jets 
are being studied. The stability characteristics are then used to compute a 
typical sound radiation pattern from a form of Lighthill's integrals. Both 
axisymmetric and non-axi symmetric disturbances are studied. Some typical 
results for axisymmetric jets are given in the attached figure, which shows 
the growth of flow instabilities in the acoustical frequency range. 


* Cambridge Hydrodynamics, Inc. 

** NASA-Langley Research Center, ANRD 

*** Kentron International, Inc., Hampton Technical Center 
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ON THE STABILITY OF ROTATING DISK FLOW 

* 

Mujeeb R. Malik 

Linear stability analysis of flow due to a rotating disk results in a 
sixth order differential equation. However, if the Coriolis and non-parallel 
effects are neglected, then we obtain the usual Orr-Sommerfel d equation. Cal- 
culations were made using this relatively simple fourth order system. For the 
purposes of transition correlation, the envelope method was employed in which 
the growth rate is maximized with respect to wave number vector and the growth 
direction is taken to be the one defined by real group velocity vector. The 
results indicate that the transition N factor is ^ 23.5 (at R = 520). The pre- 
dicted number of stationary vortices at this Reynolds number is 33 which is in 
close agreement with the experimental evidence. ^ 

The critical Reynolds number obtained by using Orr-Sommerfel d equation is 
171. However, at such a small Reynolds number it cannot be justified to 
neglect the Coriolis force and non-parallel terms. So the sixth order system 
is solved by use of Chebyshev spectral technique employing complex QR- 
algorithm. The preliminary results indicate that the critical Reynolds number 
is much less than 171 and that the growth rates are larger than those given by 
Orr-Sommerfel d equation. 


Gregory, N., Stuart, J. T., and Walker, W. S., "On the Stability of Three- 
Dimensional Boundary Layers with Application to the Flow Due to a Rotating 
Disk," Phil. Tran. Roy. Soc., A-248, 1955, pp. 155-199. 
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STABILITY OF FLOWS OVER BODIES WITH 
SUCTION THROUGH POROUS STRIPS 

All H. Nayfeh* 

One of the approaches proposed for laminar flow control is suction 
through porous strips. To determine the effectiveness of such an approach 
and to optimize the number, spacing, and flow rate through such strips, one 
needs to determine the stability of the flow over a body with suction strips. 
The first step in such an approach is the calculation of the mean flow. 

Nayfeh and El-Hady' used a nonsimilar boundary-layer code. They found a 
strip-suction distribution that is as effective as uniform suction for one 
flow condition with an increase of less than 7 percent in the overall suction. 
These results could be optimized even further. However, the computations are 
very expensive because nonsimilar boundary-layer calculations require a large 
computation time. Moreover, they fail to account for the upstream influence 
(deviation in pressure gradient) and they are limited by the magnitude of 
suction through a strip. To overcome these difficulties, a linearized triple- 
deck formulation has been developed, which is more efficient and requires 
much less computer time. Consequently, it is ideal for optimization of the 
suction distribution for laminar flow control. Since the linearized triple- 
deck equations are not expected to be valid for large suction levels and very 
wide strips, their validity limits are being determined by comparison of their 
results with those of the interacting boundary- layer equations. The resulting 
mean-flow subroutine is being coupled with the nonparallel stability codes to 
produce a design tool that can optimize the strip-suction distribution for 
laminar flow control. 


^Nayfeh, A. H., and El-Hady, N. M.: "An Evaluation of Suction Through Porous 
Strips for Laminar Flow Control." AIAA Paper 79-1494, 1979. 
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9. Atmospheric Modeling 


Methods for simulating atmospheric circulation, weather, storms. 



PRIMITIVE EQUATION SPECTRAL MODEL 
OF THE ATMOSPHERIC GENERAL CIRCULATION 

Dr. William L. Grose,* W. Thomas Blackshear,* and Richard E. Turner* 


A three-dimensional model of the atmospheric general circulation has 
been adapted to the STAR computer system. The model is based upon the 
primitive equations with dependent variables represented as truncated 
expansions of surface spherical harmonics. Vertical derivatives are ex- 
pressed as finite differences. 

In its current form, the model utilizes truncation with maximum zonal 
wave-number of 21. A realistic representation of the Earth orography has 
been implemented in the model. Thermal forcing effects due to differential 
land-sea heating have been included by an appropriate parameterization. A 
parameterization for surface drag has been adopted to prevent unrealistic 
acceleration of the surface winds. Incorporation of a scale-dependent 
dissipation was included to prevent non-linear cascade of energy to the 
small-scale waves. A time-filter was introduced to guarantee stability during 
long-term (~ 1 year) integrations. 

Current modifications are directed to incorporating diabatic heating in 
the dynamics model. Long-wave radiative processes will be treated by a 
Newtonian formulation. Ozone heating in the ultraviolet and visible 
spectrum will be included in the model stratosphere. 

Several extended runs (~ 60 days) have been made with the model to 
study the characteristics of the model and to examine the effects of the 
various parameterizations included for physical processes. Initial studies 
will be conducted with the model to study the formation of large-scale 
blocking anticyclones. 


*AESD, 304-91-00-03, 804-827-2537 
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SUMMARY OF A Z0NALLY-AVERA6ED 
CIRCULATION MODEL OF THE ATMOSPHERE 

Richard E. Turner* 


The zonally-averaged primitive equations of atmospheric circulation are 
integrated numerically to forecast the seasonal, zonally-averaged values of 
temperature, zonal wind, and meridional wind. The purpose of these computa- 
tions is to provide transport simulation for stratospheric trace constituents 
related to environmental pollution problems. 

The tendency equations for temperature, zonal wind, and meridional wind 
are formulated in the flux conserving form. Spatial derivatives are approx- 
imated by second-order accurate finite differences. Time integration is 
performed with the second-order accurate "leap frog" time differencing 
scheme. 

The turbulent transports of synoptic scale cyclones are provided by two 
techniques. The turbulent transport of zonal wind is specified from 
atmospheric data and with seasonal variation. The turbulent transport of 
quasi -conserved quantities, such as potential temperature, are computed from 
eddy diffusivity coefficients. 

Diabatic heating rates are computed from a detailed heating routine that 
accounts for atmospheric ozone, carbon dioxide, and water vapor as functions 
of season. 

In a typical experiment, the model atmosphere is started from a given 
set of initial conditions for temperature, zonal wind, and meridional wind. 

The integration is performed over several yearly heating cycles until sub- 
sequent years repeat. The resulting flow fields are used to simulate transport 
for chemically reacting atmospheric constituents. 


*AESD, 147-30-01, 804-827-2537 
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10. Aeroacoustic Methods 


Methods designed primarily for the study of aeroacoustics phenomena. 
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EXTENSION OF FARASSAT THEORY TO INCLUDE HELICOPTER BLADE VORTEX 

INTERACTION NOISE 

Danny R . Hoad* 

Impulsive noise, vrtien present, is typically rated as the most objection- 
able form of helicopter noise. Two forms of impulsive noise are known to be 
generated by most helicopters: high-speed thickness noise at the advancing blad< 
tip and low-speed partial-power, descending-flight blade-vortex interaction 
(BVI) noise. The latter form of impulsive noise is the result of the interac- 
tion of a rotor blade with the shed tip vortex from a previous blade passage. 

A theory for the determination of the acoustic pressure-time history for 
the helicopter has been developed by Farassat (ref. 1). It is based on an 
equation derived by Ffowcs Williams and Hawkings (ref. 2). The formulation 
includes a so-called thickness noise term associated with the local normal 
velocity of the moving blade and a so-called loading noise term associated with 
the local stress acting on the medium at the blade surface. The second term 
requires detailed input aerodynamic pressures on the blade. In fact, the 
sensitivity and accuracy of the formulation due to pressure variations are as 
yet not well defined. 

Figures 1 and 2 are typical pressure-time histories of helicopter noise 
without and with BVI, impulsive noise, respectively. Since a BVI primarily 
results in a local flow angle (loading) change over a short period of time, it 
can be modeled theoretically by local pressure Impulses on the blade. Using a 
simple momentum theory strip analysis for blade-span azimuth pressure calcula- 
tions and coupling the rotor geometry analysis developed by Landgrebe (ref. 3), 
the present effort is designed to provide a program to predict noise generated 
by helicopter BVI. 

One uncertainty in the calculations is the prediction of blade-vortex 
encounter. Figure 3 indicates the complexity of the rotor wake in the presence 
of a four-bladed rotor system. This figure is that predicted for a rotor 
condition at which wind-tunnel noise data are available. Using the noise data, 
it was possible to determine the precise location of the BVI, thus the 
sensitive parameters in the wake geometry program can be evaluated to determine 
their effect on the exact prediction of occurrence. With this in hand, the 
wake geometry program can be coupled with the noise prediction program to 
estimate the BVI noise characteristics as seen in figure 2. 

References: 

1. Farassat, F.: Theory of Noise Generation from Moving Bodies with an 

Application to Helicopter Rotors. NASA TR R-451, 1975. 

2. Ffowcs Williams, J. E. ; and Hawkings, D. L. : Sound Generation by Turbulence 

and Surfaces in Arbitrary Motion. Phil. Trans. Roy. Soc. London Ser. A, 

vol. 264, May 8, 1969, pp. 321-342 

3. Landgrebe, A. J. ; and Egolf, T. A.: Prediction of Helicopter Induced 

Flow Velocities Using the Rotorcraft Wake Analysis. Proceedings of 32nd 

Annual Am. Hel. Soc. Forum, A.H.S., May 1976 
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Figure 1.- Sample acoustic pressure-time history without 
blade-vortex interaction impulsive noise. 
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Figure 2.- Sample acoustic pressure-time history with 
blade-vortex interaction impulsive noise. 



Figure 3.- Sample rotor-wake geometry calculations. 
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NONLINEAR SONIC-BOOM PROPAGATION INCLUDING 
ASYMMETRIC EFFECTS AND SHOCK COALESCENCE 


Christine M. Darden* 

Experimental sonic-boom studies in recent years have been plagued with 
the necessity of keeping model size to a maximum of five or six inches. This 
constraint on model size is necessary because all available propagation methods 
assume the input to be axi symmetric, and at distances of three body lengths, 
all disturbances are assumed to be approaching their mid-field, axisymmetric 
shape. In order to acquire data at three body lengths in currently used super- 
sonic wind tunnels, the size of the model would have to be small. 

Recent results in low-boom design studies indicate great promise if design 
techniques such as wing camber, twist, etc., can be used to meet both low-boom 
area constraints and aerodynamic performance criteria. Larger models are needed 
to efficiently incorporate all of these design techniques, resulting in one model 
which can be tested for both sonic boom and aerodynamic performance. Presently 
two models are required. 

In addition to the size problem caused by available propagation methods, 
there are also problems associated with high altitudes and high Mach numbers 
(fig. 1). These methods are all based on modified linear theory and at high 
Mach numbers and high altitudes the omitted terms in the governing equations 
become significant, thus, introducing inaccuracies in the results. 

A sonic-boom propagation computer program developed under NASA Grant 
NGL-33-016-119 promises help in several of these weak areas. It will relax 
the size constraint on the models by allowing three-dimensional effects to be 
included in the input as illustrated in figure 2. A detailed description of the 
flow field in the near-field is required. Such a description is provided by 
wind-tunnel tests or by computer programs which analyze the complete flow field 
about a body. Shock and flow patterns are then extrapolated to the ground by 
means of a Modified Method of Characteristics (MMOC). This modification changes 
outgoing characteristics to a curved shape consistent with Whitham Theory and 
thus allows step sizes to increase to the size of a body length in the far field. 
Calculations are restricted to the plane of symmetry beneath the body and all 

asymmetric effects are felt through second derivatives in such as V. .. 

W 

Additional benefits of the MMOC program are listed in figure 3. Because 
of the higher order and entropy terms included in the governing equations, this 
program also provides more accurate sonic-boom predictions from high altitudes 
and high Mach numbers. One weakness of the program is its inability to calcuate 
the coalescence of embedded shocks. This prevents solutions from being 
obtained when such calculations are necessary. 

An in-house effort is now in progress to modify the program, enabling it 
to calculate shock coalescence for axisymmetric bodies and then for nonaxi- 


*HSAD, 533-01-43-04, 804-827-3181 
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symmetric bodies (fig. 4). This feature is needed in the program so that 
results from realistic low-boom transport designs can be propagated to the 
ground. Work on the axi symmetric coalescence is complete. 


• STRENGTHS 

★ MODIFIED MOC (MMOC) 

★ REAL ATMOSPHERE 

★ HIGHER-ORDER TERMS 

i( ENTROPY VARIANCE ACROSS SHOCKS 

★ NON-AXISYMMETRIC TERMS 


• WEAKNESS 

★ NO SHOCK COALESCENCE ROUTINE 


MID-FIELD DATA 



h//-3 

AXISYMMETRIC 


NEAR-FIELD DATA 



NON-AXISYMMETRIC 


Figure 1. - Features of Available 
Sonic-Boom Propaga- 
tion Methods. 


Figure 2. - Input Data Assump- 
tions. 


• AXISYMMETRIC INPUT 

— MID-FIELO INPUT DATA 
— ^ SMALL WIND TUNNEL MODELS 

— INABILITY TO MAKE PROPER USE OF CAMBER, 
TWIST AND OTHER DESIGN TOOLS. 

• MODIFIED LINEAR THEORY 

INACCURACIES AT HIGH MACH NUMBERS 
— ^ INACCURACIES AT HIGH ATTITUDES 


Figure 3. - MMOC Sonic-Boom 
Program. 




(NON-AXISYMMaRIC) 


Figure 4. - Features Being Added 
In-House. 
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VORTEX MODELING OF TURBULENT FLOWS 


Jay C. Hardin* 


A recently completed in-house study compared the flow field predicted by a 
two-dimensional discrete vortex model with experimental data for high Reynolds 
number turbulent flow over a cavity, both on a time-average and a spectral 
basis. The data were obtained during an earlier in-house assessment of sound 
generation by flow over cavities in aircraft surfaces which was conducted in 
the Aircraft Noise Reduction Laboratory. The present comparative study was 
motivated by the increasing use of vortex models to investigate the large scale 
structure in turbulent flows as well as their utility in computation of aero- 
acoustic noise generation which had been shown in previous in-house and out-of- 
house analyses. In spite of this widespread interest, only the most elementary 
comparisons between flow parameters predicted by the vortex models and actual 
data had been attempted. 

Typical results of the comparison are shown on the accompanying montage. 

In the upper left, a schematic of the experimental geometry is shown. 
Comparisons are made with flow data taken along the line marked "Measurement 
Location" which spans the intense shear layer over the cavity near its back 
wall. In the upper right, the measured mean velocity in the streamwise 
direction along this line is compared with that predicted by the vortex model. 
Note the good agreement except down in the cavity where three-dimensional 
effects become important. In the lower left, the streamwise turbulent 
intensity profile is similarly compared. Note that the levels are generally 
reasonable above the cavity but are too high by a factor of two below the top 
of the cavity. This is thought to be due to the fact that the model has no 
third dimension into which turbulent energy may be transferred. This would 
also account for the peak of the profile occurring somewhat lower in the cavity 
in the model prediction. Finally, in the lower right, the narrowband 
unnormalized streamwise velocity spectal densities are compared for the point 
Y/D = 1.35 where the measured and predicted intensities (integral of the 
spectrum) are nearly equal. Note that remarkably good agreement exists both in 
level and frequency range, with the model somewhat overpredicting at lower 
frequencies and underpredicting at high frequencies. This can be explained due 
to the lack of an energy cascade mechanism in the model. The narrowband spikes 
which appear in the spectra are the well-known self-sustained oscillation of 
flow over cavities. The model reproduces this phenomenon slightly shifted in 
frequency due to the assumption of incompressibility. 

These results are most encouraging for the use of vortex models to calcu- 
late noise and to simulate turbulent flow fields. Apparently, even in turbu- 
lent flow field, part of the flow is deterministic, "determined" by the initial 
and boundary conditions on the flow. The most crucial areas for further 
research are indicated to be the inclusion of three-dimensionality and energy 
transfer in a computationally efficient manner. 


*ANRD, 505-06-23, 804-827-2617 


240 



241 



VORTEX MODEL OF CAVITY FLOWS 


MEAN VELOCiTY PROFILE 


TURBULENT INTENSITY PROFILE 


OMPARISON OF VELOCITY SPECTRA 
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THE INTERACTION BETWEEN A SOUND PULSE AND A JET SHEAR LAYER 


I - FAR FIELD 

Lucio Maestrello*, Alvin Bayliss**, Eli Turkel 


A code has been developed to compute the fluctuating field of an acoustic 
source embedded in a jet by solving the full, time dependent Euler equations, 
linearized about a realistic spreading jet. The code is fully vectorized, and 
realistic simulation can be obtained efficiently with CDC Cyber 203. Parallel 
experimental investigations have been conducted, and these confirm the qualita- 
tive feature predicted by the code, pointing out a new feature that can be 
studied. 

An amplification of the total acoustic power output is observed when a 
source is located within the potential flow core of the jet. This amplification 
occurs in the range of frequenci.es where the local instability waves have the 
strongest growth rate. The acoustic power amplification exhibits a peak which 
is similar to that which is observed both experimentally and analytically for 
instability waves in an unexcited jet. This is particularly true when the ampli- 
fication rate is plotted as a function of Strouhal number (fz/U) based on the 
distance of the source from the nozzle, figure 1. 

These results indicate that modification of the stability characteristic 
of the jet flow can cause significant attenuation of sound. In this regard, it 
is intended to use the code to study the effect of a plug flow nozzle and 
minimize the noise radiated from a given source. 

References: 

1. Maestrello. L.; Bayliss, A.; and Turkel, E. "Experimental and Numerical 

Results on a Shear Layer Excited by a Sound Pulse" NASA TM 80183, 1979. 

2. Maestrello, L.; Bayliss, A.; and Turkel, E. "On the Interaction Between a 

Sound Pulse with the Shear Layer of an Axi symmetric Jet" Sent for publi- 
cation to the Journal of Sound and Vibration. 
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THE INTERACTION BETWEEN A SOUND PULSE AND A JET SHEAR LAYER 

II - INFLOW FIELD 

Alvin Bayliss*, Lucio Maestrello** 


The code developed in Part I has been applied to study the near and inflow 
fields for the fluctuating perturbation. This work is also being combined with 
experimental investigation. 

The numerical simulation predicts a large amplification of the inflow fluc- 
tuating field in a jet excited by a pulse. This is accompanied by an increase 
in the far-field acoustic intensity. Both phenomena occur primarily in the 
low to medium frequency range and have been observed experimentally. The fre- 
quency of the peak far-field power amplification is found to correspond to the 
frequency of the maximum growth rate of the instability wave at the position of 
the source. 

The numerical results presented here were obtained by solving a hyperbolic 
initial boundary value problem. A family of boundary conditions are presented, 
which enable the fluctuating field to be computed in a computational domain 
which is localized in the vicinity of the flow. This suggests that the inflow 
data, computed in a relatively small region near the flow, can be input as the 
shear noise sources into existing theories of aerodynamic noise to compute the 
far-field power amplification due to the shear interaction, see for example the 
variation of the Reynolds stresses and vorticity in figure 1 and the evolution 
of the inflow pressure in figure 2, 

Reference: 

1. Bayliss, A.; Maestrello, L. "The Near Field Interaction Between a Sound 
Pulse and a Jet Shear Layer" Presented at the AIAA 6th Aeroacoustics 
Conference, Hartford, Connecticut, June 1980. 
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11. Grid Generation 


Methods for constructing computational grids for solution of flow problems, 
and studies of the effects of grid on the numerical solutions. 
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DEVELOPMENT OF MESH GENERATION TECHNIQUES IN COMPUTATIONAL 
FLUID DYNAMICS FOR APPLICATIONS TO FLOW FIELD ANALYSES OF 
AERONAUTICAL COMPONENTS AND SYSTEMS WITH COMPLEX 
GEOMETRIC CONFIGURATIONS 
Peter R. Eiseman* 

In an approach to the mesh generation problem, numerous specifications of 
mesh properties became possible with the development of the multi-surface 
method [1]. With the multi -surface method, intermediate control surfaces are 
introduced between bounding surfaces so that the mesh properties can be natu- 
rally prescribed by parametric alignment. In particular, the control surfaces, 
which should not be mistaken for coordinate surfaces, define a discrete vector 
field which is interpolated, integrated, and normalized to obtain the desired 
transformation. In the basic study, only polynomial interpolates were consid- 
ered. The resulting transformations provided all of the desired boundary 
specifications for the form of the mesh. 

In the present study, the multi-surface transformation is derived with 
interpolants each of which vanishes identically outside of a small region 
around its corresponding point of interpolation. The size of each region and 
the shape of each interpolation function is chosen so that certain important 
properties such as uniformity are admissible. These local interpolants then 
lead directly to precise local controls on the form of the mesh both along the 
boundaries and internally. Although local controls are available with methods 
based upon Poisson equations and with the polynomial version of the multi-sur- 
face transformation, the capability for precise internal specifications are not 
As an example of the desired degree of precision, suppose that we wish to trans 
form a region with curved boundaries in such a way that most of the interior is 
covered by a uniform Cartesian mesh. A schematic illustration of such a 
smoothly embedded Cartesian region is given in the figure below. The smooth- 
ness of the embedding means that derivative continuity can be specified up to a 
desired level so that, for example, finite difference approximations which ex- 
tend over both Cartesian and curvilinear regions are well represented. A dis- 
crete representation of the fluid dynamic equations over the illustrated mesh 
would then smoothly change from a fully curvilinear form to a simpler Cartesian 
form over most of the mesh area. With the fully curvilinear portions con- 
strained to regions near boundaries, the simplification in problem formulation 
is clear. In this example, the previous methods of coordinate generation would 
not be very successful. However, the multi -surface transformation with local 
interpolants would be both successful and efficient in the generation of such a 
mesh. 

Reference: 

1. Eiseman, Peter R. "A Multi-Surface Method of Coordinate Generation" J. 
Computational Physics, Vol. 33, (1979), pp. 118-150. 
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Figure 1. A Coordinate System with a Rectilinear 
Embedded Cartesian System. 
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APPLICATION OF A NUMERICAL ORTHOGONAL 


COORDINATE GENERATOR TO AXISYMMETRIC BLUNT BODIES 
Randolph A. Graves, Jr.* 


In recent years, a great deal of effort has been expended to 
develop coordinate transformations and/or mesh generators for 
varying degrees of geometric complexity. These techniques range 
from the very simple algebraic manipulators to the more mathe- 
matically elegant complex variable transformations. For flow 
over blunt axisymmetric bodies, it is imperative that the boundary 
conditions on the body be represented as accurately as possible 
and in most applications, the body surface is generally used as 
one of the coordinate lines. This is called the body-oriented 
coordinate system and has found wide application, however, a 
major problem with this method is that any discontinuity in slope 
on the body makes it impossible to describe the complete flow 
field. Thus, the body oriented system cannot handle a body with 
a concavity. 

To overcome these difficulties, a simple numerical 
orthogonal coordinate generator is applied to blunt bodies both 
for total flow fields and for forebody regions only. In this 
technique, a significant improvement in handling boundary 
conditions is gained, for not only is the body a coordinate line, 
but the forebody shock for supersonic flow also is a coordinate 
line. Thus, the mesh is orthogonal on both boundaries. 

The coordinate generator has been applied to a range of 
blunt bodies both with equal and unequal coordinate spacing. 
Additionally, the origin of the outer boundary can be displaced 
to produce a "natural" compression of mesh where needed. Figure 1 
shows a typical blunt entry body with a displaced outer boundary 
and unequal spacing in the normal direction. This technique 
handles a wide range of body shapes easily and is computationally 
efficient. 

References : 

1. Graves, Jr., R. A.: Application of a Numerical Orthogonal 

Coordinate Generator to Axisymmetric Blunt Bodies, 

NASA TM 80131, Oct. 1979. 


*SSD, 506-51-13, 202-755-3277 


250 




Figure 1 - Orthogonal coordinate system about planetary entry body with 
displaced outer boundary 


251 


Algebraic Techniques for Grid Generation 


Robert E. Smith 


Boundary-fitted coordinate system are well suited for the 
application of finite difference techniques to the governing 
equations of fluid flow. Since grid lines lie on boundary 
surfaces, boundary conditions can be accurately applied and 
corresponding computer program logic minimized. An algebraic 
technique for boundary-fitted coordinate systems is one where 
an explicit functional relation establishes the transformation 
between a uniform rectangular computational domain and the 
physical domain. This approach has two advantages - rapid, 
precise, grid control and availability of analytical derivative's. 

An approach for establishing an algebraic relationship 
between a computational domain and a physical domain is based on 
defining two disconnected boundaries and a connecting function. 
This is call the "two-boundary technique" and is described in 
reference 1. The suggested connecting functions are simple 
linear and cubic parametric pol 5 momials with parameters based 
on position and derivatives at the disconnected boundaries. 

The boundaries can be defined analytically or approximately. 

A tension spline approximation to the boundaries is one 
approach. The density of the grid in specified regions is 
controlled by embedding control functions on the independent 
parametric variables for the surface definition or on the 
independent variable for the connecting function. Grids for 
spike-nosed bodies and wedge -cylinder corners described in 
Reference 2 and 3 have been generated using this technique. 

Figure 1 shows the computational domain and the physical domain 
for one plane of a 6° wedge -cylinder corner. The density solution 
obtained from a Wavier- Strokes solution is also shown. 


REFERENCES 
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Boundary-Fitted Coordinate Systems for Fluid Flow Simulation," 
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Figure 1 - Density Distribution 
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EFFECTS OF COORDINATE SYSTEMS ON THE 
NUMERICAL SOLUTION OF FLUID FLOW 
Joe F. Thompson and C. Wayne Mastin* 


The efficient solution of a fluid flow problem requires a variable grid 
spacing which depends on both accuracy and resolution requirements. The concen- 
tration of grid points may vary greatly throughout the physical region. Rapid 
change in coordinate line spacing increases the truncation error in difference 
approximations and, in the case of viscous flow problems, introduces a numeri- 
cal viscosity which may swamp the physical viscosity in high Reynolds number 
calculations. 

Taylor series expansions about grid points in the physical region can be 
used to analyse the truncation error resulting from the classical difference 
formulations as well as to develop higher order schemes on arbitrary physical 
regions. Further study could lead to developing optimal coordinate systems. 

An example of the error caused by varying the coordinate line spacing 
and the improvement in using a higher order scheme can be seen in the solution 
of the following problem. The potential function for ideal flow about a 
circular cylinder was computed numerically on the coordinate system in Figure 
1. In the traditional discretization of Laplace's equation, the physical 
region is transformed to a rectangular region and the chain rule is used to 
approximate the second order derivatives. This result is compared in Figure 2 
with an approximation obtained by a higher order method using a Taylor series 
expansion at grid points in the physical region. 

Other sources of error when using variable grids are under investigation. 
For example, there is the problem of error caused by extremely nonorthogonal 
coordinate lines. In the discretization of the transformed equation on a 
rectangular coordinate system, the coefficients may be computed numerically or 
analytically. Preliminary results suggest that the analytic computation is 
more sensitive to variations in grid spacing. When a Neumann boundary condi- 
tion must be imposed, as in the above example, there are presently two avail- 
able options. Either a one-sided difference approximation can be used or a 
false boundary can be added so that central difference approximations are 
feasible. Both methods have performed equally well in the above example and 
in the computation of ideal flow about a Joukowski airfoil. 


*Mississippi State University, Mississippi State, MS 39762, 601-325-3623 
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1. Coordinate system with a five 
jump in coordinate line spacing. 
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TIME DEPENDENT COORDINATE SYSTEMS 


FOR THE NUMERICAL SOLUTION OF FLUID FLOW 

Liviu Lustman 
Old Dominion University 

Powerful algorithms have been used recently to map complicated geometrical 
regions over standard squares or cubes. In the corresponding new coordinates, 
uniform grids may be used efficiently, as well as techniques necessitating 
simple geometry, such as spectral methods. However, the grid generation is 
usually effected only once, before starting the solution of the differential 
system. The coordinates thus obtained are appropriate for boundary layer resolu- 
tion or sharp corners, but are essentially static and probably not applicable 
to time dependent problems. It is a serious challenge to design a time dependeni 
grid, which will map a uniform computational mesh, at various times, onto sur- 
faces in the physical domain fitted to wave phenomena - sharp gradients or shock 
discontinuities propagating in space. Proper bunching of these curvilinear 
coordinates may solve the shock resolution problem, since - as numerical exper- 
ience shows - Gibbs phenomena and shock smearing usually occurs not on a fixed 
length scale, but on a fixed number of mesh points. Clearly, such an ambitious 
project must overcome several difficulties: 

1) Grid generation is time consuming, if done by the elliptic system method of 
Thompson et al , Simpler geometrical methods are hard to match to flow 
properties in the physical domain. 

2) After being expressed in the new coordinates, the transformed equations will 
contain new terms due to coordinate evolution, posing stability and accuracy 
problems, 

However, recent results of K. Miller and R. N. Miller show that a very simple 
moving node algorithm - which is moreover independent of the physical system 
to be solved - is quite efficient in one-dimensional problems. 

I propose to study similarly conceived algorithms for sample two-dimensiona' 
computations, in order to derive a gradient sensitive, time dependent coordinate 
generator, and to test its usefulness for various numerical procedures. 


References: 
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12. Additional Topics 
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DEVELOPMENT OF A RAREFIED FLOW-FIELD ANALYSIS USING 


THE DIRECT SIMULATION MONTE-CARLO METHOD 
James N. Moss* 


The in-house development of a numerical capability to analyze the flow 
about planetary and advanced space transportation vehicles for rarefied and 
transition flow conditions is in progress. Substantial impetus for this 
development comes from the recent (ref. 1) advocacy for using aerodynamic 
forces, often at low-density conditions, to capture a vehicle into a closed 
orbit and/or to modify an orbit. The two techniques being studied for future 
application are aerobraking and aerocapture. 

On the basis of previous investigations of analytical and numerical 
techniques (refs. 2 and 3), the Direct Simulation Monte-Carlo Technique of 
Bird (ref. 3) has been selected for application to the aforementioned 
problem. This method is based oh first principles in kinetic theory and 
satisfies all the conservation laws applicable to gas dynamics, surface and 
boundary interactions, and molecular collisions. It is a probabilistic 
technique which produces results that are the numerical equivalent of a 
solution for the time-dependent Boltzmann equation. The Direct Simulation 
Method is applicable to three-dimensional steady or unsteady flow with 
internal and external boundary conditions. Assumptions concerning a starting 
distribution function are not required and initial simplifying physical 
assumptions are not required. 

Figure 1 presents results obtained by Hauser and Brock (ref. 2) with 
the Direct Simulation Method for the molecular density surrounding the shuttle 
configuration. These results are for the shuttle midplane and for the free- 
stream conditions denoted in figure 1. The normalized density values are those 
for the free-stream molecules that have been reflected from the shuttle surface 
(type 2 molecules). 

References : 
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n= 8.33x10 cm ,u-7790 m sec 
» JS -2 

\*386 m , 1 / 5 , = 10 cm sec 

T*855 K . T.= 300 K 


X (meters) 


Fig, 1 The density distribution of reflected freestream molecules 
(type 2) in the Shuttle midplane, normalized by the undisturbed free 
stream density n=8.33 x 10^ cm“^. The freestream velocity 
u=7790 m sec“^, the temperature T=855 K and the mean free path 
X=386 m. The Shuttle surface temperature Ts=300 K, and the surface 
outgassing rate Vq„= 10^^ cm““ sec“l. 
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